913#ifdef SCOREP_USER_ENABLE
914#include "scorep/SCOREP_User.inc"
939 REAL,
DIMENSION(2) :: ta
942 INTEGER(mpi) :: iopnmp
953 INTEGER(mpi) :: nsecnd
954 INTEGER(mpi) :: ntsec
956 CHARACTER (LEN=24) :: chdate
957 CHARACTER (LEN=24) :: chost
959 CHARACTER (LEN=6) :: c6
960 INTEGER major, minor, patch
984 CALL getenv(
'HOSTNAME',chost)
985 IF (chost(1:1) ==
' ')
CALL getenv(
'HOST',chost)
986 WRITE(*,*)
'($Id: dd0c569a1aafb6f6eb2f26b9b9537a685639ca25 $)'
991 CALL ilaver( major,minor, patch )
992 WRITE(*,110) lapack64, major,minor, patch
993110
FORMAT(
' using LAPACK64 with ',(a),
', version ',i0,
'.',i0,
'.',i0)
995 WRITE(*,*)
'using Intel oneMKL PARDISO'
999 WRITE(*,111) __gnuc__ , __gnuc_minor__ , __gnuc_patchlevel__
1000111
FORMAT(
' compiled with gcc ',i0,
'.',i0,
'.',i0)
1003 WRITE(*,111) __pgic__ , __pgic_minor__ , __pgic_patchlevel__
1004111
FORMAT(
' compiled with pgi ',i0,
'.',i0,
'.',i0)
1006#ifdef SCOREP_USER_ENABLE
1007 WRITE(*,*)
'instrumenting Score-P user regions'
1010 WRITE(*,*)
' < Millepede II-P starting ... ',chdate
1011 WRITE(*,*)
' ',chost
1014 WRITE(8,*)
'($Id: dd0c569a1aafb6f6eb2f26b9b9537a685639ca25 $)'
1016 WRITE(8,*)
'Log-file Millepede II-P ', chdate
1017 WRITE(8,*)
' ', chost
1019 CALL peend(-1,
'Still running or crashed')
1028 WRITE(*,*)
'!!! Checking input only, no calculation of a solution !!!'
1029 WRITE(8,*)
'!!! Checking input only, no calculation of a solution !!!'
1048 CALL getenv(
'OMP_NUM_THREADS',c6)
1050 CALL getenv(lapack64//
'_NUM_THREADS',c6)
1052 IF (c6(1:1) ==
' ')
THEN
1054 WRITE(*,*)
'Number of threads for LAPACK: unkown (empty OMP_NUM_THREADS)'
1056 WRITE(*,*)
'Number of threads for LAPACK: unkown (empty ',lapack64//
'_NUM_THREADS)'
1059 WRITE(*,*)
'Number of threads for LAPACK: ', c6
1079 CALL mvopen(lun,
'millepede.his')
1085 CALL mvopen(1,
'mpdebug.txt')
1089 times(0)=rstext-rstp
1095 times(1)=rloop1-rstext
1099 WRITE(8,*)
'Chi square cut equiv 3 st.dev applied ...'
1100 WRITE(8,*)
' in first iteration with factor',
chicut
1101 WRITE(8,*)
' in second iteration with factor',
chirem
1102 WRITE(8,*)
' (reduced by sqrt in next iterations)'
1106 WRITE(8,*)
'Down-weighting of outliers in',
lhuber,
' iterations'
1107 WRITE(8,*)
'Cut on downweight fraction',
dwcut
1111 times(2)=rloop2-rloop1
1116 CALL peend(5,
'Ended without solution (empty constraints)')
1118 CALL peend(0,
'Ended normally')
1155 CALL gmpdef(6,1,
'log10(#records) vs file number')
1156 CALL gmpdef(7,1,
'final rejection fraction vs file number')
1158 'final <Chi^2/Ndf> from accepted local fits vs file number')
1159 CALL gmpdef(9,1,
'<Ndf> from accepted local fits vs file number')
1165 rej=real(nrc-
jfd(kfl),mps)/real(nrc,mps)
1166 CALL gmpxy(6,real(kfl,mps),log10(real(nrc,mps)))
1167 CALL gmpxy(7,real(kfl,mps),rej)
1169 IF (
jfd(kfl) > 0)
THEN
1170 c2ndf=
cfd(kfl)/real(
jfd(kfl),mps)
1171 CALL gmpxy(8,real(kfl,mps),c2ndf)
1172 andf=real(
dfd(kfl),mps)/real(
jfd(kfl),mps)
1173 CALL gmpxy(9,real(kfl,mps),andf)
1190 WRITE(*,*)
'Misalignment test wire chamber'
1193 CALL hmpdef( 9,-0.0015,+0.0015,
'True - fitted displacement')
1194 CALL hmpdef(10,-0.0015,+0.0015,
'True - fitted Vdrift')
1200 sums(1)=sums(1)+diff
1201 sums(2)=sums(2)+diff*diff
1203 sums(3)=sums(3)+diff
1204 sums(4)=sums(4)+diff*diff
1206 sums(1)=0.01_mpd*sums(1)
1207 sums(2)=sqrt(0.01_mpd*sums(2))
1208 sums(3)=0.01_mpd*sums(3)
1209 sums(4)=sqrt(0.01_mpd*sums(4))
1210 WRITE(*,143)
'Parameters 1 - 100: mean =',sums(1),
'rms =',sums(2)
1211 WRITE(*,143)
'Parameters 101 - 200: mean =',sums(3),
'rms =',sums(4)
1212143
FORMAT(6x,a28,f9.6,3x,a5,f9.6)
1215 WRITE(*,*)
' I label simulated fitted diff'
1216 WRITE(*,*)
' -------------------------------------------- '
1236 WRITE(*,*)
'Misalignment test Si tracker'
1239 CALL hmpdef( 9,-0.0025,+0.0025,
'True - fitted displacement X')
1240 CALL hmpdef(10,-0.025,+0.025,
'True - fitted displacement Y')
1251 sums(1)=sums(1)+1.0_mpd
1252 sums(2)=sums(2)+diff
1253 sums(3)=sums(3)+diff*diff
1258 err=sqrt(abs(gmati))
1260 sums(7)=sums(7)+1.0_mpd
1261 sums(8)=sums(8)+diff
1262 sums(9)=sums(9)+diff*diff
1265 IF (mod(i,3) == 1)
THEN
1269 sums(4)=sums(4)+1.0_mpd
1270 sums(5)=sums(5)+diff
1271 sums(6)=sums(6)+diff*diff
1276 err=sqrt(abs(gmati))
1278 sums(7)=sums(7)+1.0_mpd
1279 sums(8)=sums(8)+diff
1280 sums(9)=sums(9)+diff*diff
1285 sums(2)=sums(2)/sums(1)
1286 sums(3)=sqrt(sums(3)/sums(1))
1287 sums(5)=sums(5)/sums(4)
1288 sums(6)=sqrt(sums(6)/sums(4))
1289 WRITE(*,143)
'Parameters 1 - 500: mean =',sums(2),
'rms =',sums(3)
1290 WRITE(*,143)
'Parameters 501 - 700: mean =',sums(5),
'rms =',sums(6)
1291 IF (sums(7) > 0.5_mpd)
THEN
1292 sums(8)=sums(8)/sums(7)
1293 sums(9)=sqrt(sums(9)/sums(7))
1294 WRITE(*,143)
'Parameter pulls, all: mean =',sums(8),
'rms =',sums(9)
1298 WRITE(*,*)
' I label simulated fitted diff'
1299 WRITE(*,*)
' -------------------------------------------- '
1311 IF (mod(i,3) == 1)
THEN
1331 WRITE(8,*)
'Record',
nrec1,
' has largest residual:',
value1
1334 WRITE(8,*)
'Record',
nrec2,
' has largest Chi^2/Ndf:',
value2
1338 WRITE(8,*)
'Record',
nrec3,
' is first with error (rank deficit/NaN)'
1342 WRITE(8,*)
'In total 3 +',
nloopn,
' loops through the data files'
1344 WRITE(8,*)
'In total 2 +',
nloopn,
' loops through the data files'
1348 WRITE(8,*)
'In total ',
mnrsit,
' internal MINRES iterations'
1356 ntsec=nint(deltat,mpi)
1357 CALL sechms(deltat,nhour,minut,secnd)
1358 nsecnd=nint(secnd,mpi)
1359 WRITE(8,*)
'Total time =',ntsec,
' seconds =',nhour,
' h',minut, &
1360 ' m',nsecnd,
' seconds'
1362 WRITE(8,*)
'end ', chdate
1369 IF(sum(
nrejec) /= 0)
THEN
1371 WRITE(8,*)
'Data records rejected in last iteration: '
1378 WRITE(*,*)
' < Millepede II-P ending ... ', chdate
1385106
FORMAT(
' PARDISO dyn. memory allocation: ',f11.6,
' GB')
1395 WRITE(*,*)
'Postprocessing:'
1405102
FORMAT(2x,i4,i10,2x,3f10.5)
1406103
FORMAT(
' Times [in sec] for text processing',f12.3/ &
1409 ' func. value ',f12.3,
' *',f4.0/ &
1410 ' func. value, global matrix, solution',f12.3,
' *',f4.0/ &
1411 ' new solution',f12.3,
' *',f4.0/)
1412105
FORMAT(
' Peak dynamic memory allocation: ',f11.6,
' GB')
1432 INTEGER(mpi) :: istop
1433 INTEGER(mpi) :: itgbi
1434 INTEGER(mpi) :: itgbl
1436 INTEGER(mpi) :: itnlim
1437 INTEGER(mpi) :: nout
1439 INTEGER(mpi),
INTENT(IN) :: ivgbi
1476 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1480 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1483 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1489 err=sqrt(abs(real(gmati,mps)))
1490 IF(gmati < 0.0_mpd) err=-err
1491 diag=matij(ivgbi,ivgbi)
1492 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps)
1494101
FORMAT(1x,
' label parameter presigma differ', &
1495 ' Error gcor^2 iit'/ 1x,
'---------',2x,5(
'-----------'),2x,
'----')
1496102
FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1516 INTEGER(mpi) :: istop
1517 INTEGER(mpi) :: itgbi
1518 INTEGER(mpi) :: itgbl
1520 INTEGER(mpi) :: itnlim
1521 INTEGER(mpi) :: nout
1523 INTEGER(mpi),
INTENT(IN) :: ivgbi
1552 mxxnrm = real(
nagb,mpd)/sqrt(epsilon(mxxnrm))
1554 trcond = 1.0_mpd/epsilon(trcond)
1555 ELSE IF(
mrmode == 2)
THEN
1563 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1567 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1571 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1578 err=sqrt(abs(real(gmati,mps)))
1579 IF(gmati < 0.0_mpd) err=-err
1580 diag=matij(ivgbi,ivgbi)
1581 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps)
1583101
FORMAT(1x,
' label parameter presigma differ', &
1584 ' Error gcor^2 iit'/ 1x,
'---------',2x,5(
'-----------'),2x,
'----')
1585102
FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1598 INTEGER(mpi) :: icgb
1599 INTEGER(mpi) :: irhs
1600 INTEGER(mpi) :: itgbi
1601 INTEGER(mpi) :: ivgb
1603 INTEGER(mpi) :: jcgb
1605 INTEGER(mpi) :: label
1607 INTEGER(mpi) :: inone
1610 REAL(mpd) :: drhs(4)
1611 INTEGER(mpi) :: idrh (4)
1636 IF(abs(rhs) > climit)
THEN
1642 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1651 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1654 WRITE(*,102)
' Constraints: only equation values >', climit,
' are printed'
1655101
FORMAT(
' ',4(i6,g11.3))
1656102
FORMAT(a,g11.2,a)
1669 INTEGER(mpi) :: icgb
1670 INTEGER(mpi) :: icgrp
1671 INTEGER(mpi) :: ioff
1672 INTEGER(mpi) :: itgbi
1674 INTEGER(mpi) :: jcgb
1675 INTEGER(mpi) :: label
1676 INTEGER(mpi) :: labelf
1677 INTEGER(mpi) :: labell
1678 INTEGER(mpi) :: last
1679 INTEGER(mpi) :: line1
1680 INTEGER(mpi) :: ncon
1681 INTEGER(mpi) :: ndiff
1682 INTEGER(mpi) :: npar
1683 INTEGER(mpi) :: inone
1684 INTEGER(mpi) :: itype
1685 INTEGER(mpi) :: ncgbd
1686 INTEGER(mpi) :: ncgbr
1687 INTEGER(mpi) :: ncgbw
1688 INTEGER(mpi) :: ncgrpd
1689 INTEGER(mpi) :: ncgrpr
1690 INTEGER(mpi) :: next
1692 INTEGER(mpl):: length
1693 INTEGER(mpl) :: rows
1695 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecParConsOffsets
1696 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecParConsList
1697 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsParOffsets
1698 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsParList
1699 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: matConsGroupIndex
1712 IF(last < 0.AND.label < 0)
THEN
1715 IF(itype == 2) ncgbw=ncgbw+1
1722 IF(label > 0.AND.itype == 2)
THEN
1729 IF (ncgbw == 0)
THEN
1730 WRITE(*,*)
'GRPCON:',
ncgb,
' constraints found in steering files'
1732 WRITE(*,*)
'GRPCON:',
ncgb,
' constraints found in steering files,',ncgbw,
'weighted'
1737 length=
ncgb+1; rows=3
1750 CALL mpalloc(matconsgroupindex,rows,length,
'group index for constraint (I)')
1753 CALL mpalloc(vecconsparoffsets,length,
'offsets for global par list for cons. (I)')
1755 CALL mpalloc(vecparconsoffsets,length,
'offsets for cons. list for global par. (I)')
1756 vecparconsoffsets(1)=0
1758 vecparconsoffsets(i+1)=vecparconsoffsets(i)+
globalparcons(i)
1762 length=vecparconsoffsets(
ntgb+1)
1763 CALL mpalloc(vecconsparlist,length,
'global par. list for constraint (I)')
1764 CALL mpalloc(vecparconslist,length,
'constraint list for global par. (I)')
1769 vecconsparoffsets(1)=ioff
1781 vecparconslist(vecparconsoffsets(itgbi)+
globalparcons(itgbi))=icgb
1783 vecconsparlist(ioff+npar)=itgbi
1789 CALL sort1k(vecconsparlist(ioff+1),npar)
1793 next=vecconsparlist(ioff+j)
1794 IF (next /= last)
THEN
1796 vecconsparlist(ioff+ndiff) = next
1805 vecconsparoffsets(icgb+1)=ioff
1818 print *,
' Constraint #parameters first par. last par. first line'
1826 npar=vecconsparoffsets(icgb+1)-vecconsparoffsets(icgb)
1830 print *, jcgb, npar, labelf, labell, line1
1833 icgrp=matconsgroupindex(1,icgb)
1834 IF (icgrp == 0)
THEN
1836 DO i=vecconsparoffsets(icgb)+1, vecconsparoffsets(icgb+1)
1837 itgbi=vecconsparlist(i)
1839 DO j=vecparconsoffsets(itgbi)+1,vecparconsoffsets(itgbi+1)
1840 icgrp=matconsgroupindex(1,vecparconslist(j))
1846 IF (icgrp == 0)
THEN
1853 matconsgroupindex(2,icgb)=jcgb
1854 matconsgroupindex(3,icgb)=icgb
1855 DO i=vecconsparoffsets(icgb)+1, vecconsparoffsets(icgb+1)
1856 itgbi=vecconsparlist(i)
1859 DO j=vecparconsoffsets(itgbi)+1,vecparconsoffsets(itgbi+1)
1860 matconsgroupindex(1,vecparconslist(j))=icgrp
1864 WRITE(*,*)
'GRPCON:',
ncgrp,
' disjoint constraints groups built'
1872 icgb=matconsgroupindex(3,jcgb)
1877 icgrp=matconsgroupindex(1,jcgb)
1897 print *,
' cons.group first con. first par. last par. #cons #par'
1908 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell, ncon, npar
1911 IF (ncon == npar)
THEN
1918 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell,
' : cons.group resolved'
1933 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell,
' : cons.group redundant'
1938 IF (ncgrpr > 0)
THEN
1939 WRITE(*,*)
'GRPCON:',ncgbr,
' redundancy constraints in ', ncgrpr,
' groups resolved'
1943 IF (ncgrpd > 0)
THEN
1944 WRITE(*,*)
'GRPCON:',ncgbd,
' redundancy constraints in ', ncgrpd,
' groups detected'
1967 INTEGER(mpi) :: icgb
1968 INTEGER(mpi) :: icgrp
1969 INTEGER(mpi) :: ifrst
1970 INTEGER(mpi) :: ilast
1971 INTEGER(mpi) :: isblck
1972 INTEGER(mpi) :: itgbi
1973 INTEGER(mpi) :: ivgb
1975 INTEGER(mpi) :: jcgb
1976 INTEGER(mpi) :: jfrst
1977 INTEGER(mpi) :: label
1978 INTEGER(mpi) :: labelf
1979 INTEGER(mpi) :: labell
1980 INTEGER(mpi) :: ncon
1981 INTEGER(mpi) :: ngrp
1982 INTEGER(mpi) :: npar
1983 INTEGER(mpi) :: ncnmxb
1984 INTEGER(mpi) :: ncnmxg
1985 INTEGER(mpi) :: nprmxb
1986 INTEGER(mpi) :: nprmxg
1987 INTEGER(mpi) :: inone
1988 INTEGER(mpi) :: nvar
1990 INTEGER(mpl):: length
1991 INTEGER(mpl) :: rows
1993 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: matConsGroupIndex
2006 length=
ncgrp+1; rows=3
2011 CALL mpalloc(matconsgroupindex,rows,length,
'group index for constraint (I)')
2048 IF (nvar > 0 .OR.
iskpec == 0)
THEN
2051 matconsgroupindex(1,
ncgb)=ngrp+1
2052 matconsgroupindex(2,
ncgb)=icgb
2053 matconsgroupindex(3,
ncgb)=nvar
2056 IF (
ncgb > ncon) ngrp=ngrp+1
2062 WRITE(*,*)
'PRPCON:',
ncgbe,
' empty constraints skipped'
2064 WRITE(*,*)
'PRPCON:',
ncgbe,
' empty constraints detected, to be fixed !!!'
2065 WRITE(*,*)
' (use option "skipemptycons" to skip those)'
2070 WRITE(*,*)
'!!! Switch to "-C" (checking input only), no calculation of a solution !!!'
2071 WRITE(8,*)
'!!! Switch to "-C" (checking input only), no calculation of a solution !!!'
2076 WRITE(*,*)
'PRPCON:',
ncgb,
' constraints accepted'
2079 IF(
ncgb == 0)
RETURN
2086 icgb=matconsgroupindex(2,jcgb)
2091 icgrp=matconsgroupindex(1,jcgb)
2117 WRITE(*,*)
' Cons. sorted index #var.par. first line first label last label'
2118 WRITE(*,*)
' Cons. group index first cons. last cons. first label last label'
2119 WRITE(*,*)
' Cons. block index first group last group first label last label'
2125 nvar=matconsgroupindex(3,jcgb)
2129 WRITE(*,*)
' Cons. sorted', jcgb, nvar, &
2132 WRITE(*,*)
' Cons. sorted', jcgb,
' empty (0)', &
2141 WRITE(*,*)
' Cons. group ', icgrp,
matconsgroups(1,icgrp), &
2154 DO i=icgrp,isblck,-1
2168 WRITE(*,*)
' Cons. block ',
ncblck, isblck, icgrp, labelf, labell
2186 ncnmxb=max(ncnmxb,ncon)
2187 nprmxb=max(nprmxb,npar)
2212 ncnmxg=max(ncnmxg,ncon)
2213 nprmxg=max(nprmxg,npar)
2248 WRITE(*,*)
'PRPCON: constraints split into ',
ncgrp,
'(disjoint) groups,'
2249 WRITE(*,*)
' groups combined into ',
ncblck,
'(non overlapping) blocks'
2250 WRITE(*,*)
' max group size (cons., par.) ', ncnmxg, nprmxg
2251 WRITE(*,*)
' max block size (cons., par.) ', ncnmxb, nprmxb
2269 INTEGER(mpi) :: icgb
2270 INTEGER(mpi) :: icgrp
2272 INTEGER(mpi) :: ifirst
2273 INTEGER(mpi) :: ilast
2274 INTEGER(mpl) :: ioffc
2275 INTEGER(mpl) :: ioffp
2276 INTEGER(mpi) :: irank
2277 INTEGER(mpi) :: ipar0
2278 INTEGER(mpi) :: itgbi
2279 INTEGER(mpi) :: ivgb
2281 INTEGER(mpi) :: jcgb
2283 INTEGER(mpi) :: label
2284 INTEGER(mpi) :: ncon
2285 INTEGER(mpi) :: npar
2286 INTEGER(mpi) :: nrank
2287 INTEGER(mpi) :: inone
2292 INTEGER(mpl):: length
2293 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: matConstraintsT
2294 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: auxVectorD
2295 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: auxVectorI
2299 IF(
ncgb == 0)
RETURN
2308 CALL mpalloc(auxvectori,length,
'auxiliary array (I)')
2309 CALL mpalloc(auxvectord,length,
'auxiliary array (D)')
2312 CALL mpalloc(matconstraintst,length,
'transposed matrix of constraints (blocks)')
2313 matconstraintst=0.0_mpd
2326 WRITE(*,*)
' Constraint group, #con, rank', icgrp, ncon, 0,
' (empty)'
2329 DO jcgb=ifirst,ilast
2341 IF(ivgb > 0) matconstraintst(int(jcgb-ifirst,mpl)*int(npar,mpl)+ivgb-ipar0+ioffc)= &
2342 matconstraintst(int(jcgb-ifirst,mpl)*int(npar,mpl)+ivgb-ipar0+ioffc)+factr
2349 DO ll=ioffc+1,ioffc+npar
2355 matconstraintst(int(i-1,mpl)*int(npar,mpl)+ll)* &
2356 matconstraintst(int(j-1,mpl)*int(npar,mpl)+ll)
2362 IF (
icheck > 1 .OR. irank < ncon)
THEN
2363 WRITE(*,*)
' Constraint group, #con, rank', icgrp, ncon, irank
2364 IF (irank < ncon)
THEN
2365 WRITE(*,*)
' .. rank deficit !! '
2366 WRITE(*,*)
' E.g. fix all parameters and remove all constraints related to label ', &
2371 ioffc=ioffc+int(npar,mpl)*int(ncon,mpl)
2378 WRITE(*,*)
'Rank of product matrix of constraints is',nrank, &
2379 ' for',
ncgb,
' constraint equations'
2380 WRITE(8,*)
'Rank of product matrix of constraints is',nrank, &
2381 ' for',
ncgb,
' constraint equations'
2382 IF(nrank <
ncgb)
THEN
2383 WRITE(*,*)
'Warning: insufficient constraint equations!'
2384 WRITE(8,*)
'Warning: insufficient constraint equations!'
2387 WRITE(*,*)
' --> enforcing SUBITO mode'
2388 WRITE(8,*)
' --> enforcing SUBITO mode'
2395 print *,
'QL decomposition of constraints matrix'
2398 WRITE(
lunlog,*)
'QL decomposition of constraints matrix'
2410 CALL lpqldec(matconstraintst,evmin,evmax)
2414 print *,
' largest |eigenvalue| of L: ', evmax
2415 print *,
' smallest |eigenvalue| of L: ', evmin
2416 IF (evmin == 0.0_mpd.AND.
icheck == 0)
THEN
2417 CALL peend(27,
'Aborted, singular QL decomposition of constraints matrix')
2418 stop
'FEASMA: stopping due to singular QL decomposition of constraints matrix'
2444 INTEGER(mpi) :: icgb
2445 INTEGER(mpi) :: icgrp
2446 INTEGER(mpi) :: iter
2447 INTEGER(mpi) :: itgbi
2448 INTEGER(mpi) :: ivgb
2449 INTEGER(mpi) :: ieblck
2450 INTEGER(mpi) :: isblck
2451 INTEGER(mpi) :: ifirst
2452 INTEGER(mpi) :: ilast
2454 INTEGER(mpi) :: jcgb
2455 INTEGER(mpi) :: label
2456 INTEGER(mpi) :: inone
2457 INTEGER(mpi) :: ncon
2459 REAL(mps),
INTENT(IN) :: concut
2460 INTEGER(mpi),
INTENT(OUT) :: iact
2467 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: vecCorrections
2471 IF(
ncgb == 0)
RETURN
2501 sum1=sqrt(sum1/real(
ncgb,mpd))
2502 sum2=sum2/real(
ncgb,mpd)
2504 IF(iter == 1.AND.sum1 < concut)
RETURN
2506 IF(iter == 1.AND.
ncgb <= 12)
THEN
2508 WRITE(*,*)
'Constraint equation discrepancies:'
2510101
FORMAT(4x,4(i5,g12.4))
2512103
FORMAT(10x,
' Cut on rms value is',g8.1)
2517 WRITE(*,*)
'Improve constraints'
2521 WRITE(*,102) iter,sum1,sum2,sum3
2522102
FORMAT(i6,
' rms',g12.4,
' avrg_abs',g12.4,
' max_abs',g12.4)
2524 CALL mpalloc(veccorrections,int(
nvgb,mpl),
'constraint corrections')
2525 veccorrections=0.0_mpd
2533 ieblck=isblck+(ncon*(ncon+1))/2
2602 INTEGER(mpi) :: iact
2603 INTEGER(mpi) :: ierrc
2604 INTEGER(mpi) :: ierrf
2605 INTEGER(mpi) :: ioffp
2607 INTEGER(mpi) :: ithr
2608 INTEGER(mpi) :: jfile
2609 INTEGER(mpi) :: jrec
2611 INTEGER(mpi) :: kfile
2614 INTEGER(mpi) :: mpri
2616 INTEGER(mpi) :: nact
2617 INTEGER(mpi) :: nbuf
2618 INTEGER(mpi) :: ndata
2619 INTEGER(mpi) :: noff
2620 INTEGER(mpi) :: noffs
2621 INTEGER(mpi) :: npointer
2622 INTEGER(mpi) :: npri
2626 INTEGER(mpi) :: nrpr
2627 INTEGER(mpi) :: nthr
2628 INTEGER(mpi) :: ntot
2629 INTEGER(mpi) :: maxRecordSize
2630 INTEGER(mpi) :: maxRecordFile
2632 INTEGER(mpi),
INTENT(OUT) :: more
2642 CHARACTER (LEN=7) :: cfile
2647 SUBROUTINE readc(bufferD, bufferF, bufferI, bufferLength, lun, err)
BIND(c)
2649 REAL(c_double),
DIMENSION(*),
INTENT(OUT) :: bufferD
2650 REAL(c_float),
DIMENSION(*),
INTENT(OUT) :: bufferF
2651 INTEGER(c_int),
DIMENSION(*),
INTENT(OUT) :: bufferI
2652 INTEGER(c_int),
INTENT(INOUT) :: bufferLength
2653 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
2654 INTEGER(c_int),
INTENT(OUT) :: err
2655 END SUBROUTINE readc
2661 DATA npri / 0 /, mpri / 1000 /
2709 files:
DO WHILE (jfile > 0)
2713 CALL binopn(kfile,ithr,ios)
2719 IF(kfile <=
nfilf)
THEN
2742 eof=(ierrc <= 0.AND.ierrc /= -4)
2743 IF(eof.AND.ierrc < 0)
THEN
2744 WRITE(*,*)
'Read error for binary Cfile', kfile,
'record', jrec+1,
':', ierrc
2745 WRITE(8,*)
'Read error for binary Cfile', kfile,
'record', jrec+1,
':', ierrc
2747 WRITE(cfile,
'(I7)') kfile
2748 CALL peend(18,
'Aborted, read error(s) for binary file ' // cfile)
2749 stop
'PEREAD: stopping due to read errors (bad record, wrong file type?)'
2751 IF (
kfd(1,jfile) == 1)
THEN
2757 IF(eof)
EXIT records
2762 xfd(jfile)=max(
xfd(jfile),n)
2779 IF ((noff+nr+2+
ndimbuf >= ndata*(iact+1)).OR.(nbuf >= npointer))
EXIT files
2794 IF (
kfd(1,jfile) == 1)
THEN
2795 print *,
'PEREAD: file ', kfile,
'read the first time, found',jrec,
' records'
2799 IF (-
kfd(1,jfile) /= jrec)
THEN
2800 WRITE(cfile,
'(I7)') kfile
2801 CALL peend(19,
'Aborted, binary file modified (length) ' // cfile)
2802 stop
'PEREAD: file modified (length)'
2853 IF (
kfd(1,jfile) == 1)
kfd(1,jfile)=-nrc
2872 DO WHILE (
nloopn == 0.AND.ntot >= nrpr)
2873 WRITE(*,*)
' Record ',nrpr
2874 IF (nrpr < 100000)
THEN
2883 IF (npri == 1)
WRITE(*,100)
2885100
FORMAT(/
' PeRead records active file' &
2886 /
' total block threads number')
2887101
FORMAT(
' PeRead',4i10)
2900 IF (
xfd(k) > maxrecordsize)
THEN
2901 maxrecordsize=
xfd(k)
2904 dw=real(-
kfd(1,k),mpd)
2907 ds1=ds1+dw*real(
wfd(k),mpd)
2908 ds2=ds2+dw*real(
wfd(k)**2,mpd)
2910 print *,
'PEREAD: file ', maxrecordfile,
'with max record size ', maxrecordsize
2911 IF (
nfilw > 0.AND.ds0 > 0.0_mpd)
THEN
2915 WRITE(lun,177)
nfilw,real(ds1,mps),real(ds2,mps)
2916177
FORMAT(/
' !!!!!',i4,
' weighted binary files', &
2917 /
' !!!!! mean, variance of weights =',2g12.4)
2935179
FORMAT(/
' Read cache usage (#blocks, #records, ', &
2936 'min,max records/block'/17x,i10,i12,2i10)
2954 INTEGER(mpi),
INTENT(IN) :: mode
2956 INTEGER(mpi) :: ibuf
2957 INTEGER(mpi) :: ichunk
2959 INTEGER(mpi) :: itgbi
2965 INTEGER(mpi),
PARAMETER :: maxbad = 100
2966 INTEGER(mpi) :: nbad
2967 INTEGER(mpi) :: nerr
2968 INTEGER(mpi) :: inone
2986 CALL isjajb(nst,ist,ja,jb,jsp)
3005#ifdef SCOREP_USER_ENABLE
3006 scorep_user_region_by_name_begin(
"UR_peprep", scorep_user_region_type_common)
3011 CALL pechk(ibuf,nerr)
3014 IF(nbad >= maxbad)
EXIT
3019 CALL isjajb(nst,ist,ja,jb,jsp)
3032 CALL peend(20,
'Aborted, bad binary records')
3033 stop
'PEREAD: stopping due to bad records'
3036#ifdef SCOREP_USER_ENABLE
3037 scorep_user_region_by_name_end(
"UR_peprep")
3057 INTEGER(mpi) :: ioff
3064 INTEGER(mpi),
INTENT(IN) :: ibuf
3065 INTEGER(mpi),
INTENT(OUT) :: nerr
3074 outer:
DO WHILE(is < nst)
3079 IF(is > nst)
EXIT outer
3085 IF(is > nst)
EXIT outer
3102100
FORMAT(
' PEREAD: record ', i8,
' in file ',i6,
' is broken !!!')
3112101
FORMAT(
' PEREAD: record ', i8,
' in file ',i6,
' contains ', i6,
' NaNs !!!')
3128 INTEGER(mpi) :: ibuf
3129 INTEGER(mpi) :: ichunk
3130 INTEGER(mpi) :: iproc
3131 INTEGER(mpi) :: ioff
3132 INTEGER(mpi) :: ioffbi
3134 INTEGER(mpi) :: itgbi
3139 INTEGER(mpi) :: nalg
3140 INTEGER(mpi) :: neqna
3143 INTEGER(mpi) :: nzero
3144 INTEGER(mpi) :: inone
3145 INTEGER(mpl) :: length
3180 CALL isjajb(nst,ist,ja,jb,jsp)
3211 CALL isjajb(nst,ist,ja,jb,jsp)
3239#ifdef SCOREP_USER_ENABLE
3240 scorep_user_region_by_name_begin(
"UR_pepgrp", scorep_user_region_type_common)
3249 CALL pargrp(ist+1,ist+nnz)
3261#ifdef SCOREP_USER_ENABLE
3262 scorep_user_region_by_name_end(
"UR_pepgrp")
3281 INTEGER(mpi) :: istart
3282 INTEGER(mpi) :: itgbi
3284 INTEGER(mpi) :: jstart
3285 INTEGER(mpi) :: jtgbi
3286 INTEGER(mpi) :: lstart
3287 INTEGER(mpi) :: ltgbi
3289 INTEGER(mpi),
INTENT(IN) :: inds
3290 INTEGER(mpi),
INTENT(IN) :: inde
3292 IF (inds > inde)
RETURN
3301 IF (istart == 0)
THEN
3302 IF (itgbi /= ltgbi+1)
THEN
3305 IF (lstart == 0)
THEN
3320 IF (istart /= jstart)
THEN
3331 IF (itgbi /= ltgbi+1)
THEN
3346 IF (istart /= jstart)
THEN
3398 INTEGER(mpi),
INTENT(IN) :: nst
3399 INTEGER(mpi),
INTENT(IN OUT) :: is
3400 INTEGER(mpi),
INTENT(OUT) :: ja
3401 INTEGER(mpi),
INTENT(OUT) :: jb
3402 INTEGER(mpi),
INTENT(OUT) :: jsp
3410 IF(is >= nst)
RETURN
3464 INTEGER(mpi) :: ioffb
3466 INTEGER(mpi) :: itgbi
3467 INTEGER(mpi) :: itgbij
3468 INTEGER(mpi) :: itgbik
3469 INTEGER(mpi) :: ivgb
3470 INTEGER(mpi) :: ivgbij
3471 INTEGER(mpi) :: ivgbik
3474 INTEGER(mpi) :: lastit
3476 INTEGER(mpi) :: ncrit
3477 INTEGER(mpi) :: ngras
3478 INTEGER(mpi) :: nparl
3480 INTEGER(mpl) :: nrej
3481 INTEGER(mpi) :: inone
3482 INTEGER(mpi) :: ilow
3483 INTEGER(mpi) :: nlow
3484 INTEGER(mpi) :: nzero
3504 CALL gmpdef(1,4,
'Function value in iterations')
3506 CALL gmpdef(2,3,
'Number of MINRES steps vs iteration nr')
3508 CALL hmpdef( 5,0.0,0.0,
'Number of degrees of freedom')
3509 CALL hmpdef(11,0.0,0.0,
'Number of local parameters')
3510 CALL hmpdef(16,0.0,24.0,
'LOG10(cond(band part decomp.)) local fit ')
3511 CALL hmpdef(23,0.0,0.0,
'SQRT of diagonal elements without presigma')
3512 CALL hmpdef(24,0.0,0.0,
'Log10 of off-diagonal elements')
3513 CALL hmpdef(25,0.0,0.0,
'Relative individual pre-sigma')
3514 CALL hmpdef(26,0.0,0.0,
'Relative global pre-sigma')
3519 'Normalized residuals of single (global) measurement')
3521 'Normalized residuals of single (local) measurement')
3523 'Pulls of single (global) measurement')
3525 'Pulls of single (local) measurement')
3526 CALL hmpdef(4,0.0,0.0,
'Chi^2/Ndf after local fit')
3527 CALL gmpdef(4,5,
'location, dispersion (res.) vs record nr')
3528 CALL gmpdef(5,5,
'location, dispersion (pull) vs record nr')
3542 IF(
nloopn == 2)
CALL hmpdef(6,0.0,0.0,
'Down-weight fraction')
3545 IF(
iterat /= lastit)
THEN
3553 ELSE IF(
iterat >= 1)
THEN
3601 WRITE(*,111) nparl,ncrit,usei,used,peaki,peakd
3602111
FORMAT(
' Write cache usage (#flush,#overrun,<levels>,', &
3603 'peak(levels))'/2i7,
',',4(f6.1,
'%'))
3611 ELSE IF (
mbandw > 0)
THEN
3643 print *,
" ... warning ..."
3644 print *,
" global parameters with too few (< MREQENA) accepted entries: ", nlow
3648 IF(
icalcm == 1 .AND. nzero > 0)
THEN
3650 WRITE(*,*)
'Warning: the rank defect of the symmetric',
nfgb, &
3651 '-by-',
nfgb,
' matrix is ',
ndefec,
' (should be zero).'
3652 WRITE(lun,*)
'Warning: the rank defect of the symmetric',
nfgb, &
3653 '-by-',
nfgb,
' matrix is ',
ndefec,
' (should be zero).'
3656 WRITE(*,*)
' --> enforcing SUBITO mode'
3657 WRITE(lun,*)
' --> enforcing SUBITO mode'
3667 elmt=real(matij(i,i),mps)
3668 IF(elmt > 0.0)
CALL hmpent(23,1.0/sqrt(elmt))
3702 CALL addsums(1, adder, 0, 1.0_mpl)
3716 IF(sgm > 0.0) weight=1.0/sgm**2
3732 IF(itgbij /= 0)
THEN
3738 IF (factrj == 0.0_mpd) cycle
3748 IF(
icalcm == 1.AND.ivgbij > 0)
THEN
3755 IF(ivgbij > 0.AND.ivgbik > 0)
THEN
3756 CALL mupdat(ivgbij,ivgbik,weight*factrj*factrk)
3762 adder=weight*dsum**2
3763 CALL addsums(1, adder, 1, 1.0_mpl)
3781 ratae=max(-99.9,ratae)
3790 WRITE(*,*)
'Data records rejected in initial loop:'
3835 WRITE(lun,*)
' === local fits have bordered band matrix structure ==='
3836 IF (
nbndr(1) > 0 )
WRITE(lun,101)
' NBNDR',
nbndr(1),
'number of records (upper/left border)'
3837 IF (
nbndr(2) > 0 )
WRITE(lun,101)
' NBNDR',
nbndr(2),
'number of records (lower/right border)'
3838 WRITE(lun,101)
' NBDRX',
nbdrx,
'max border size'
3839 WRITE(lun,101)
' NBNDX',
nbndx,
'max band width'
3848101
FORMAT(1x,a8,
' =',i14,
' = ',a)
3863 INTEGER(mpi),
INTENT(IN) :: lunp
3868101
FORMAT(
' it fc',
' fcn_value dfcn_exp slpr costh iit st', &
3869 ' ls step cutf',1x,
'rejects hhmmss FMS')
3870102
FORMAT(
' -- --',
' ----------- -------- ---- ----- --- --', &
3871 ' -- ----- ----',1x,
'------- ------ ---')
3887 INTEGER(mpl) :: nrej
3888 INTEGER(mpi) :: nsecnd
3892 REAL(mps) :: slopes(3)
3893 REAL(mps) :: steps(3)
3894 REAL,
DIMENSION(2) :: ta
3897 INTEGER(mpi),
INTENT(IN) :: lunp
3899 CHARACTER (LEN=4):: ccalcm(4)
3900 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3904 IF(nrej > 9999999) nrej=9999999
3908 nsecnd=nint(secnd,mpi)
3917 CALL ptlopt(nfa,ma,slopes,steps)
3918 ratae=max(-99.9,min(99.9,slopes(2)/slopes(1)))
3924103
FORMAT(i3,i3,e12.5,38x,f5.1, 1x,i7, i3,i2.2,i2.2,a4)
3925104
FORMAT(i3,i3,e12.5,1x,e8.2,f6.3,f6.3,i5,2i3,f6.3,f5.1, &
3926 1x,i7, i3,i2.2,i2.2,a4)
3927105
FORMAT(i3,i3,e12.5,1x,e8.2,12x,i5,i3,9x,f5.1, &
3928 1x,i7, i3,i2.2,i2.2,a4)
3941 INTEGER(mpi) :: minut
3943 INTEGER(mpi) :: nhour
3944 INTEGER(mpl) :: nrej
3945 INTEGER(mpi) :: nsecnd
3949 REAL(mps) :: slopes(3)
3950 REAL(mps) :: steps(3)
3951 REAL,
DIMENSION(2) :: ta
3954 INTEGER(mpi),
INTENT(IN) :: lunp
3955 CHARACTER (LEN=4):: ccalcm(4)
3956 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3960 IF(nrej > 9999999) nrej=9999999
3964 nsecnd=nint(secnd,mpi)
3968 CALL ptlopt(nfa,ma,slopes,steps)
3969 ratae=abs(slopes(2)/slopes(1))
3974104
FORMAT(3x,i3,e12.5,9x, 35x, i7, i3,i2.2,i2.2,a4)
3975105
FORMAT(3x,i3,e12.5,9x, f6.3,14x,i3,f6.3,6x, i7, i3,i2.2,i2.2,a4)
3989 INTEGER(mpi) :: nsecnd
3992 REAL,
DIMENSION(2) :: ta
3995 INTEGER(mpi),
INTENT(IN) :: lunp
3996 CHARACTER (LEN=4):: ccalcm(4)
3997 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
4002 nsecnd=nint(secnd,mpi)
4004 WRITE(lunp,106) nhour,minut,nsecnd,ccalcm(
lcalcm)
4005106
FORMAT(69x,i3,i2.2,i2.2,a4)
4015 INTEGER(mpi) :: lunit
4017 WRITE(lunit,102)
'Explanation of iteration table'
4018 WRITE(lunit,102)
'=============================='
4019 WRITE(lunit,101)
'it', &
4020 'iteration number. Global parameters are improved for it > 0.'
4021 WRITE(lunit,102)
'First function evaluation is called iteraton 0.'
4022 WRITE(lunit,101)
'fc',
'number of function evaluations.'
4023 WRITE(lunit,101)
'fcn_value',
'value of 2 x Likelihood function (LF).'
4024 WRITE(lunit,102)
'The final value is the chi^2 value of the fit and should'
4025 WRITE(lunit,102)
'be about equal to the NDF (see below).'
4026 WRITE(lunit,101)
'dfcn_exp', &
4027 'expected reduction of the value of the Likelihood function (LF)'
4028 WRITE(lunit,101)
'slpr',
'ratio of the actual slope to inital slope.'
4029 WRITE(lunit,101)
'costh', &
4030 'cosine of angle between search direction and -gradient'
4032 WRITE(lunit,101)
'iit', &
4033 'number of internal iterations in MINRES algorithm'
4034 WRITE(lunit,101)
'st',
'stop code of MINRES algorithm'
4035 WRITE(lunit,102)
'< 0: rhs is very special, with beta2 = 0'
4036 WRITE(lunit,102)
'= 0: rhs b = 0, i.e. the exact solution is x = 0'
4037 WRITE(lunit,102)
'= 1 requested accuracy achieved, as determined by rtol'
4038 WRITE(lunit,102)
'= 2 reasonable accuracy achieved, given eps'
4039 WRITE(lunit,102)
'= 3 x has converged to an eigenvector'
4040 WRITE(lunit,102)
'= 4 matrix ill-conditioned (Acond has exceeded 0.1/eps)'
4041 WRITE(lunit,102)
'= 5 the iteration limit was reached'
4042 WRITE(lunit,102)
'= 6 Matrix x vector does not define a symmetric matrix'
4043 WRITE(lunit,102)
'= 7 Preconditioner does not define a symmetric matrix'
4044 ELSEIF (
metsol == 5)
THEN
4045 WRITE(lunit,101)
'iit', &
4046 'number of internal iterations in MINRES-QLP algorithm'
4047 WRITE(lunit,101)
'st',
'stop code of MINRES-QLP algorithm'
4048 WRITE(lunit,102)
'= 1: beta_{k+1} < eps, iteration k is the final Lanczos step.'
4049 WRITE(lunit,102)
'= 2: beta2 = 0. If M = I, b and x are eigenvectors of A.'
4050 WRITE(lunit,102)
'= 3: beta1 = 0. The exact solution is x = 0.'
4051 WRITE(lunit,102)
'= 4: A solution to (poss. singular) Ax = b found, given rtol.'
4052 WRITE(lunit,102)
'= 5: A solution to (poss. singular) Ax = b found, given eps.'
4053 WRITE(lunit,102)
'= 6: Pseudoinverse solution for singular LS problem, given rtol.'
4054 WRITE(lunit,102)
'= 7: Pseudoinverse solution for singular LS problem, given eps.'
4055 WRITE(lunit,102)
'= 8: The iteration limit was reached.'
4056 WRITE(lunit,102)
'= 9: The operator defined by Aprod appears to be unsymmetric.'
4057 WRITE(lunit,102)
'=10: The operator defined by Msolve appears to be unsymmetric.'
4058 WRITE(lunit,102)
'=11: The operator defined by Msolve appears to be indefinite.'
4059 WRITE(lunit,102)
'=12: xnorm has exceeded maxxnorm or will exceed it next iteration.'
4060 WRITE(lunit,102)
'=13: Acond has exceeded Acondlim or 0.1/eps.'
4061 WRITE(lunit,102)
'=14: Least-squares problem but no converged solution yet.'
4062 WRITE(lunit,102)
'=15: A null vector obtained, given rtol.'
4064 WRITE(lunit,101)
'ls',
'line search info'
4065 WRITE(lunit,102)
'< 0 recalculate function'
4066 WRITE(lunit,102)
'= 0: N or STP lt 0 or step not descending'
4067 WRITE(lunit,102)
'= 1: Linesearch convergence conditions reached'
4068 WRITE(lunit,102)
'= 2: interval of uncertainty at lower limit'
4069 WRITE(lunit,102)
'= 3: max nr of line search calls reached'
4070 WRITE(lunit,102)
'= 4: step at the lower bound'
4071 WRITE(lunit,102)
'= 5: step at the upper bound'
4072 WRITE(lunit,102)
'= 6: rounding error limitation'
4073 WRITE(lunit,101)
'step', &
4074 'the factor for the Newton step during the line search. Usually'
4076 'a value of 1 gives a sufficient reduction of the LF. Oherwise'
4077 WRITE(lunit,102)
'other step values are tried.'
4078 WRITE(lunit,101)
'cutf', &
4079 'cut factor. Local fits are rejected, if their chi^2 value'
4081 'is larger than the 3-sigma chi^2 value times the cut factor.'
4082 WRITE(lunit,102)
'A cut factor of 1 is used finally, but initially a larger'
4083 WRITE(lunit,102)
'factor may be used. A value of 0.0 means no cut.'
4084 WRITE(lunit,101)
'rejects',
'total number of rejected local fits.'
4085 WRITE(lunit,101)
'hmmsec',
'the time in hours (h), minutes (mm) and seconds.'
4086 WRITE(lunit,101)
'FMS',
'calculation of Function value, Matrix, Solution.'
4089101
FORMAT(a9,
' = ',a)
4106 INTEGER(mpi),
INTENT(IN) :: i
4107 INTEGER(mpi),
INTENT(IN) :: j
4108 REAL(mpd),
INTENT(IN) :: add
4110 INTEGER(mpl):: ijadd
4111 INTEGER(mpl):: ijcsr3
4116 IF(i <= 0.OR.j <= 0.OR. add == 0.0_mpd)
RETURN
4137 ELSE IF(
matsto == 2)
THEN
4165 CALL peend(23,
'Aborted, bad matrix index')
4166 stop
'mupdat: bad index'
4191 INTEGER(mpi),
INTENT(IN) :: i
4192 INTEGER(mpi),
INTENT(IN) :: j1
4193 INTEGER(mpi),
INTENT(IN) :: j2
4194 INTEGER(mpi),
INTENT(IN) :: il
4195 INTEGER(mpi),
INTENT(IN) :: jl
4196 INTEGER(mpi),
INTENT(IN) :: n
4197 REAL(mpd),
INTENT(IN) :: sub((n*n+n)/2)
4204 INTEGER(mpi):: iblast
4205 INTEGER(mpi):: iblock
4211 INTEGER(mpi):: jblast
4212 INTEGER(mpi):: jblock
4222 IF(i <= 0.OR.j1 <= 0.OR.j2 > i)
RETURN
4236 print *,
' MGUPDT: ij=0', i,j1,j2,il,jl,ij,lr,iprc,
matsto
4243 ijl=(k*k-k)/2+jl+ir-ia1
4260 ijl=(k*k-k)/2+jl+ir-ia1
4264 IF (jblock /= jblast.OR.iblock /= iblast)
THEN
4288 CALL ijpgrp(i,j1,ij,lr,iprc)
4290 print *,
' MGUPDT: ij=0', i,j1,j2,il,jl,ij,lr,iprc,
matsto
4354SUBROUTINE loopbf(nrej,numfil,naccf,chi2f,ndff)
4381 INTEGER(mpi) :: ibuf
4382 INTEGER(mpi) :: ichunk
4383 INTEGER(mpl) :: icmn
4384 INTEGER(mpl) :: icost
4386 INTEGER(mpi) :: idiag
4388 INTEGER(mpi) :: iext
4396 INTEGER(mpi) :: imeas
4399 INTEGER(mpi) :: ioffb
4400 INTEGER(mpi) :: ioffc
4401 INTEGER(mpi) :: ioffd
4402 INTEGER(mpi) :: ioffe
4403 INTEGER(mpi) :: ioffi
4404 INTEGER(mpi) :: ioffq
4405 INTEGER(mpi) :: iprc
4406 INTEGER(mpi) :: iprcnx
4407 INTEGER(mpi) :: iprdbg
4408 INTEGER(mpi) :: iproc
4409 INTEGER(mpi) :: irbin
4410 INTEGER(mpi) :: isize
4412 INTEGER(mpi) :: iter
4413 INTEGER(mpi) :: itgbi
4414 INTEGER(mpi) :: ivgbj
4415 INTEGER(mpi) :: ivgbk
4416 INTEGER(mpi) :: ivpgrp
4426 INTEGER(mpi) :: joffd
4427 INTEGER(mpi) :: joffi
4428 INTEGER(mpi) :: jproc
4432 INTEGER(mpi) :: kbdr
4433 INTEGER(mpi) :: kbdrx
4434 INTEGER(mpi) :: kbnd
4437 INTEGER(mpi) :: lvpgrp
4438 INTEGER(mpi) :: mbdr
4439 INTEGER(mpi) :: mbnd
4440 INTEGER(mpi) :: mside
4441 INTEGER(mpi) :: nalc
4442 INTEGER(mpi) :: nalg
4446 INTEGER(mpi) :: ndown
4448 INTEGER(mpi) :: nfred
4449 INTEGER(mpi) :: nfrei
4451 INTEGER(mpi) :: nprdbg
4452 INTEGER(mpi) :: nrank
4455 INTEGER(mpi) :: nter
4456 INTEGER(mpi) :: nweig
4457 INTEGER(mpi) :: ngrp
4458 INTEGER(mpi) :: npar
4460 INTEGER(mpl),
INTENT(IN OUT) :: nrej(6)
4461 INTEGER(mpi),
INTENT(IN) :: numfil
4462 INTEGER(mpi),
INTENT(IN OUT) :: naccf(numfil)
4463 REAL(mps),
INTENT(IN OUT) :: chi2f(numfil)
4464 INTEGER(mpi),
INTENT(IN OUT) :: ndff(numfil)
4474 INTEGER(mpi) :: ijprec
4481 CHARACTER (LEN=3):: chast
4482 DATA chuber/1.345_mpd/
4483 DATA cauchy/2.3849_mpd/
4531 IF(
nloopn == 1.AND.mod(nrc,100000_mpl) == 0)
THEN
4532 WRITE(*,*)
'Record',nrc,
' ... still reading'
4533 IF(
monpg1>0)
WRITE(
lunlog,*)
'Record',nrc,
' ... still reading'
4543 IF(nrc ==
nrecpr) lprnt=.true.
4544 IF(nrc ==
nrecp2) lprnt=.true.
4545 IF(nrc ==
nrecer) lprnt=.true.
4550 IF (nprdbg == 1) iprdbg=iproc
4551 IF (iproc /= iprdbg) lprnt=.false.
4556 WRITE(1,*)
'------------------ Loop',
nloopn, &
4557 ': Printout for record',nrc,iproc
4568 CALL isjajb(nst,ist,ja,jb,jsp)
4570 IF(imeas == 0)
WRITE(1,1121)
45751121
FORMAT(/
'Measured value and local derivatives'/ &
4576 ' i measured std_dev index...derivative ...')
45771122
FORMAT(i3,2g12.4,3(i3,g12.4)/(27x,3(i3,g12.4)))
4583 CALL isjajb(nst,ist,ja,jb,jsp)
4585 IF(imeas == 0)
WRITE(1,1123)
45971123
FORMAT(/
'Global derivatives'/ &
4598 ' i label gindex vindex derivative ...')
45991124
FORMAT(i3,2(i9,i7,i7,g12.4)/(3x,2(i9,i7,i7,g12.4)))
46001125
FORMAT(i3,2(i9,i7,i7,g12.4))
4610 WRITE(1,*)
'Data corrections using values of global parameters'
4611 WRITE(1,*)
'=================================================='
4621 CALL isjajb(nst,ist,ja,jb,jsp)
4653101
FORMAT(
' index measvalue corrvalue sigma')
4654102
FORMAT(i6,2x,2g12.4,
' +-',g12.4)
4656 IF(nalc <= 0)
GO TO 90
4658 ngg=(nalg*nalg+nalg)/2
4671 IF (ivpgrp /= lvpgrp)
THEN
4679 IF (npar /= nalg)
THEN
4680 print *,
' mismatch of number of global parameters ', nrc, nalg, npar, ngrp
4690 CALL peend(35,
'Aborted, mismatch of number of global parameters')
4691 stop
' mismatch of number of global parameters '
4718 DO WHILE(iter < nter)
4723 WRITE(1,*)
'Outlier-suppression iteration',iter,
' of',nter
4724 WRITE(1,*)
'=========================================='
4734 DO i=1,(nalc*nalc+nalc)/2
4745 wght =1.0_mpd/rerr**2
4750 IF(abs(resid) > chuber*rerr)
THEN
4751 wght=wght*chuber*rerr/abs(resid)
4755 wght=wght/(1.0+(resid/rerr/cauchy)**2)
4759 IF(lprnt.AND.iter /= 1.AND.nter /= 1)
THEN
4761 IF(abs(resid) > chuber*rerr) chast=
'* '
4762 IF(abs(resid) > 3.0*rerr) chast=
'** '
4763 IF(abs(resid) > 6.0*rerr) chast=
'***'
4764 IF(imeas == 0)
WRITE(1,*)
'Second loop: accumulate'
4765 IF(imeas == 0)
WRITE(1,103)
4770 WRITE(1,104) imeas,rmeas,resid,rerr,r1,chast,r2
4772103
FORMAT(
' index corrvalue residuum sigma', &
4774104
FORMAT(i6,2x,2g12.4,
' +-',g12.4,f7.2,1x,a3,f8.2)
4789 im=min(nalc+1-ij,nalc+1-ik)
4796 IF (iter == 1.AND.nalc > 5.AND.
lfitbb > 0)
THEN
4799 icmn=int(nalc,mpl)**3
4805 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
4806 IF (icost < icmn)
THEN
4818 kbdr=max(kbdr,
ibandh(k+nalc))
4819 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
4820 IF (icost < icmn)
THEN
4844 IF (
icalcm == 1.OR.lprnt) inv=2
4845 IF (mside == 1)
THEN
4853 IF (evdmin > 0.0_mpl) cndl10=log10(real(evdmax/evdmin,mps))
4863 IF ((.NOT.(
blvec(k) <= 0.0_mpd)).AND. (.NOT.(
blvec(k) > 0.0_mpd))) nan=nan+1
4868 WRITE(1,*)
'Parameter determination:',nalc,
' parameters,',
' rank=',nrank
4869 WRITE(1,*)
'-----------------------'
4870 IF(ndown /= 0)
WRITE(1,*)
' ',ndown,
' data down-weighted'
4886 wght =1.0_mpd/rerr**2
4910 IF (0.999999_mpd/wght > dvar)
THEN
4911 pull=rmeas/sqrt(1.0_mpd/wght-dvar)
4914 CALL hmpent(13,real(pull,mps))
4915 CALL gmpms(5,rec,real(pull,mps))
4917 CALL hmpent(14,real(pull,mps))
4936 IF(iter == 1.AND.jb < ist.AND.lhist) &
4937 CALL gmpms(4,rec,real(rmeas/rerr,mps))
4939 dchi2=wght*rmeas*rmeas
4944 IF(abs(resid) > chuber*rerr)
THEN
4945 wght=wght*chuber*rerr/abs(resid)
4946 dchi2=2.0*chuber*(abs(resid)/rerr-0.5*chuber)
4949 wght=wght/(1.0_mpd+(resid/rerr/cauchy)**2)
4950 dchi2=log(1.0_mpd+(resid/rerr/cauchy)**2)*cauchy**2
4960 IF(abs(resid) > chuber*rerr) chast=
'* '
4961 IF(abs(resid) > 3.0*rerr) chast=
'** '
4962 IF(abs(resid) > 6.0*rerr) chast=
'***'
4963 IF(imeas == 0)
WRITE(1,*)
'Third loop: single residuals'
4964 IF(imeas == 0)
WRITE(1,105)
4968 IF(resid < 0.0) r1=-r1
4969 IF(resid < 0.0) r2=-r2
4972105
FORMAT(
' index corrvalue residuum sigma', &
4974106
FORMAT(i6,2x,2g12.4,
' +-',g12.4,f7.2,1x,a3,f8.2)
4976 IF(iter == nter)
THEN
4978 resmax=max(resmax,abs(rmeas)/rerr)
4981 IF(iter == 1.AND.lhist)
THEN
4983 CALL hmpent( 3,real(rmeas/rerr,mps))
4985 CALL hmpent(12,real(rmeas/rerr,mps))
4992 resing=(real(nweig,mps)-real(suwt,mps))/real(nweig,mps)
4994 IF(iter == 1)
CALL hmpent( 5,real(ndf,mps))
4995 IF(iter == 1)
CALL hmpent(11,real(nalc,mps))
5000 WRITE(1,*)
'Chi^2=',summ,
' at',ndf,
' degrees of freedom: ', &
5001 '3-sigma limit is',chindl(3,ndf)*real(ndf,mps)
5002 WRITE(1,*) suwt,
' is sum of factors, compared to',nweig, &
5003 ' Downweight fraction:',resing
5009 WRITE(1,*)
' NaNs ', nalc, nrank, nan
5010 WRITE(1,*)
' ---> rejected!'
5012 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-1 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
5015 IF(nrank /= nalc)
THEN
5019 WRITE(1,*)
' rank deficit', nalc, nrank
5020 WRITE(1,*)
' ---> rejected!'
5022 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-2 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
5029 WRITE(1,*)
' too large condition(band part) ', nalc, nrank, cndl10
5030 WRITE(1,*)
' ---> rejected!'
5032 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-3 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
5038 WRITE(1,*)
' Ndf<=0', nalc, nrank, ndf
5039 WRITE(1,*)
' ---> rejected!'
5041 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-4 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
5045 chndf=real(summ/real(ndf,mpd),mps)
5047 IF(iter == 1.AND.lhist)
CALL hmpent(4,chndf)
5060 chichi=chindl(3,ndf)*real(ndf,mps)
5064 IF(summ >
chhuge*chichi)
THEN
5067 WRITE(1,*)
' ---> rejected!'
5075 IF(summ > chlimt)
THEN
5077 WRITE(1,*)
' ---> rejected!'
5081 CALL addsums(iproc+1, dchi2, ndf, dw1)
5091 CALL addsums(iproc+1, dchi2, ndf, dw1)
5095 WRITE(1,*)
' ---> rejected!'
5108 naccf(kfl)=naccf(kfl)+1
5109 ndff(kfl) =ndff(kfl) +ndf
5110 chi2f(kfl)=chi2f(kfl)+chndf
5123 wght =1.0_mpd/rerr**2
5124 dchi2=wght*rmeas*rmeas
5128 IF(resid > chuber*rerr)
THEN
5129 wght=wght*chuber*rerr/resid
5130 dchi2=2.0*chuber*(resid/rerr-0.5*chuber)
5180 CALL addsums(iproc+1, summ, ndf, dw1)
5236 IF (nfred < 0.OR.nfrei < 0)
THEN
5263 iprcnx=ijprec(i,jnx)
5265 IF (.NOT.((jnx == j+1).AND.(iprc == iprcnx)))
THEN
5279 joffd=joffd+(il*il-il)/2
5291 WRITE(1,*)
'------------------ End of printout for record',nrc
5348 iprcnx=ijprec(i,jnx)
5350 IF (.NOT.((jnx == j+1).AND.(iprc == iprcnx)))
THEN
5365 joffd=joffd+(il*il-il)/2
5401 INTEGER(mpi),
INTENT(IN) :: lun
5403 IF (
nrejec(1)>0)
WRITE(lun,*)
nrejec(1),
' (local solution contains NaNs)'
5404 IF (
nrejec(2)>0)
WRITE(lun,*)
nrejec(2),
' (local matrix with rank deficit)'
5405 IF (
nrejec(3)>0)
WRITE(lun,*)
nrejec(3),
' (local matrix with ill condition)'
5406 IF (
nrejec(4)>0)
WRITE(lun,*)
nrejec(4),
' (local fit with Ndf=0)'
5407 IF (
nrejec(5)>0)
WRITE(lun,*)
nrejec(5),
' (local fit with huge Chi2(Ndf))'
5408 IF (
nrejec(6)>0)
WRITE(lun,*)
nrejec(6),
' (local fit with large Chi2(Ndf))'
5434 INTEGER(mpi) :: icom
5435 INTEGER(mpl) :: icount
5439 INTEGER(mpi) :: imin
5440 INTEGER(mpi) :: iprlim
5441 INTEGER(mpi) :: isub
5442 INTEGER(mpi) :: itgbi
5443 INTEGER(mpi) :: itgbl
5444 INTEGER(mpi) :: ivgbi
5446 INTEGER(mpi) :: label
5454 INTEGER(mpi) :: labele(3)
5455 REAL(mps):: compnt(3)
5460 CALL mvopen(lup,
'millepede.res')
5463 WRITE(*,*)
' Result of fit for global parameters'
5464 WRITE(*,*)
' ==================================='
5469 WRITE(lup,*)
'Parameter ! first 3 elements per line are', &
5470 ' significant (if used as input)'
5488 err=sqrt(abs(real(gmati,mps)))
5489 IF(gmati < 0.0_mpd) err=-err
5492 IF(gmati*diag > 0.0_mpd)
THEN
5493 gcor2=1.0_mpd-1.0_mpd/(gmati*diag)
5494 IF(gcor2 >= 0.0_mpd.AND.gcor2 <= 1.0_mpd) gcor=real(sqrt(gcor2),mps)
5499 IF(lowstat) icount=-(icount+1)
5501 IF(itgbi <= iprlim)
THEN
5515 ELSE IF(itgbi == iprlim+1)
THEN
5516 WRITE(* ,*)
'... (further printout suppressed, but see log file)'
5530 ELSE IF (
igcorr /= 0)
THEN
5548 CALL mvopen(lup,
'millepede.eve')
5558 DO isub=0,min(15,imin-1)
5580 WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
5584 IF(iev /= 0)
WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
5592101
FORMAT(1x,
' label parameter presigma differ', &
5593 ' error'/ 1x,
'-----------',4x,4(
'-------------'))
5594102
FORMAT(i10,2x,4g14.5,f8.3)
5595103
FORMAT(3(i11,f11.7,2x))
5596110
FORMAT(i10,2x,2g14.5,28x,i12)
5597111
FORMAT(i10,2x,3g14.5,14x,i12)
5598112
FORMAT(i10,2x,4g14.5,i12)
5620 INTEGER(mpi) :: icom
5621 INTEGER(mpl) :: icount
5622 INTEGER(mpi) :: ifrst
5623 INTEGER(mpi) :: ilast
5624 INTEGER(mpi) :: inext
5625 INTEGER(mpi) :: itgbi
5626 INTEGER(mpi) :: itgbl
5627 INTEGER(mpi) :: itpgrp
5628 INTEGER(mpi) :: ivgbi
5630 INTEGER(mpi) :: icgrp
5631 INTEGER(mpi) :: ipgrp
5633 INTEGER(mpi) :: jpgrp
5635 INTEGER(mpi) :: label1
5636 INTEGER(mpi) :: label2
5637 INTEGER(mpi) :: ncon
5638 INTEGER(mpi) :: npair
5639 INTEGER(mpi) :: nstep
5642 INTEGER(mpl):: length
5644 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecPairedParGroups
5647 SUBROUTINE ggbmap(ipgrp,npair,npgrp)
5649 INTEGER(mpi),
INTENT(IN) :: ipgrp
5650 INTEGER(mpi),
INTENT(OUT) :: npair
5651 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npgrp
5659 CALL mvopen(lup,
'millepede.res')
5660 WRITE(lup,*)
'*** Results of checking input only, no solution performed ***'
5661 WRITE(lup,*)
'! === global parameters ==='
5662 WRITE(lup,*)
'! fixed-1: by pre-sigma, -2: by entries cut, -3: by iterated entries cut'
5664 WRITE(lup,*)
'! Label Value Pre-sigma SkippedEntries Cons. group Status '
5666 WRITE(lup,*)
'! Label Value Pre-sigma Entries Cons. group Status '
5682 IF (ivgbi <= 0)
THEN
5684 IF (ivgbi == -4)
THEN
5685 WRITE(lup,116) c1,itgbl,par,presig,icount,icgrp
5687 WRITE(lup,110) c1,itgbl,par,presig,icount,icgrp,ivgbi
5691 WRITE(lup,111) c1,itgbl,par,presig,icount,icgrp
5697 WRITE(lup,*)
'!.Appearance statistics '
5698 WRITE(lup,*)
'!. Label First file and record Last file and record #files #paired-par'
5701 IF (itpgrp > 0)
THEN
5709 WRITE(lup,*)
'* === constraint groups ==='
5711 WRITE(lup,*)
'* Group #Cons. Entries First label Last label'
5713 WRITE(lup,*)
'* Group #Cons. Entries First label Last label Paired label range'
5715 CALL mpalloc(vecpairedpargroups,length,
'paired global parameter groups (I)')
5727 IF (
icheck > 1 .AND. label1 > 0)
THEN
5731 vecpairedpargroups(npair+1)=0
5735 jpgrp=vecpairedpargroups(j)
5739 IF (ifrst /= 0.AND.inext /= (ilast+nstep))
THEN
5742 WRITE(lup,114) label1, label2
5747 IF (ifrst == 0) ifrst=inext
5754 IF (jpgrp == vecpairedpargroups(j+1)-1)
THEN
5759 IF (ifrst /= 0)
THEN
5762 WRITE(lup,114) label1, label2
5768 WRITE(lup,*)
'*.Appearance statistics '
5769 WRITE(lup,*)
'*. Group First file and record Last file and record #files'
5779110
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' fixed',i2)
5780111
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' variable')
5781112
FORMAT(
' !.',i10,6i11)
5782113
FORMAT(
' * ',i6,i8,3i12)
5783114
FORMAT(
' *:',48x,i12,
' ..',i12)
5784115
FORMAT(
' *.',i10,5i11)
5785116
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' redundant')
5815 INTEGER(mpi) :: iproc
5825 INTEGER(mpi),
INTENT(IN) :: n
5826 INTEGER(mpl),
INTENT(IN) :: l
5827 REAL(mpd),
INTENT(IN) :: x(n)
5828 INTEGER(mpi),
INTENT(IN) :: is
5829 INTEGER(mpi),
INTENT(IN) :: ie
5830 REAL(mpd),
INTENT(OUT) :: b(n)
5835 INTEGER(mpl) :: indij
5836 INTEGER(mpl) :: indid
5838 INTEGER(mpi) :: ichunk
5874 CALL peend(24,
'Aborted, vector/matrix size mismatch')
5875 stop
'AVPRDS: mismatched vector and matrix'
5895 IF (ia2 <= ib2) b(ia2:ib2)=b(ia2:ib2)+
globalmatd(ia2:ib2)*x(ia2:ib2)
5906 IF (i <= ie.AND.i >= is)
THEN
5914 IF (j <= ie.AND.j >= is)
THEN
5930 IF (ja2 <= jb2)
THEN
5933 b(i)=b(i)+dot_product(
globalmatd(indij+lj+ja2-ja:indij+lj+jb2-ja),x(ja2:jb2))
5937 IF (
mextnd == 0.AND.ia2 <= ib2)
THEN
5940 b(j)=b(j)+dot_product(
globalmatd(indij+lj+jn*(ia2-ia):indij+lj+jn*(ib2-ia):jn),x(ia2:ib2))
5959 IF (i <= ie.AND.i >= is)
THEN
5967 IF (j <= ie.AND.j >= is)
THEN
5983 IF (ja2 <= jb2)
THEN
5986 b(i)=b(i)+dot_product(real(
globalmatf(indij+lj+ja2-ja:indij+lj+jb2-ja),mpd),x(ja2:jb2))
5990 IF (
mextnd == 0.AND.ia2 <= ib2)
THEN
5993 b(j)=b(j)+dot_product(real(
globalmatf(indij+lj+jn*(ia2-ia):indij+lj+jn*(ib2-ia):jn),mpd),x(ia2:ib2))
6027 INTEGER(mpi) :: iproc
6035 INTEGER(mpi),
INTENT(IN) :: n
6036 INTEGER(mpl),
INTENT(IN) :: l
6037 REAL(mpd),
INTENT(IN) :: x(n)
6038 REAL(mpd),
INTENT(OUT) :: b(n)
6043 INTEGER(mpl) :: indij
6044 INTEGER(mpl) :: indid
6046 INTEGER(mpi) :: ichunk
6075 CALL peend(24,
'Aborted, vector/matrix size mismatch')
6076 stop
'AVPRD0: mismatched vector and matrix'
6122 b(i)=b(i)+dot_product(
globalmatd(indij+lj:indij+lj+jn-1),x(ja:jb))
6128 b(j)=b(j)+dot_product(
globalmatd(indij+lj:indij+jn*in:jn),x(ia:ib))
6164 b(i)=b(i)+dot_product(real(
globalmatf(indij+lj:indij+lj+jn-1),mpd),x(ja:jb))
6170 b(j)=b(j)+dot_product(real(
globalmatf(indij+lj:indij+jn*in:jn),mpd),x(ia:ib))
6194 INTEGER(mpi) :: ispc
6205 INTEGER(mpl) :: indid
6206 INTEGER(mpl),
DIMENSION(12) :: icount
6213 icount(4)=huge(icount(4))
6214 icount(7)=huge(icount(7))
6215 icount(10)=huge(icount(10))
6225 ir=ipg+(ispc-1)*(
napgrp+1)
6232 icount(1)=icount(1)+in
6233 icount(2)=icount(2)+in*ku
6239 icount(3)=icount(3)+1
6240 icount(4)=min(icount(4),jn)
6241 icount(5)=icount(5)+jn
6242 icount(6)=max(icount(6),jn)
6243 icount(7)=min(icount(7),in)
6244 icount(8)=icount(8)+in
6245 icount(9)=max(icount(9),in)
6246 icount(10)=min(icount(10),in*jn)
6247 icount(11)=icount(11)+in*jn
6248 icount(12)=max(icount(12),in*jn)
6254 WRITE(*,*)
"analysis of sparsity structure"
6255 IF (icount(1) > 0)
THEN
6256 WRITE(*,101)
"rows without compression/blocks ", icount(1)
6257 WRITE(*,101)
" contained elements ", icount(2)
6259 WRITE(*,101)
"number of block matrices ", icount(3)
6260 avg=real(icount(5),mps)/real(icount(3),mps)
6261 WRITE(*,101)
"number of columns (min,mean,max) ", icount(4), avg, icount(6)
6262 avg=real(icount(8),mps)/real(icount(3),mps)
6263 WRITE(*,101)
"number of rows (min,mean,max) ", icount(7), avg, icount(9)
6264 avg=real(icount(11),mps)/real(icount(3),mps)
6265 WRITE(*,101)
"number of elements (min,mean,max) ", icount(10), avg, icount(12)
6266101
FORMAT(2x,a34,i10,f10.3,i10)
6285 INTEGER(mpi),
INTENT(IN) :: n
6286 REAL(mpd),
INTENT(IN) :: x(n)
6287 REAL(mpd),
INTENT(OUT) :: b(n)
6292 CALL peend(24,
'Aborted, vector/matrix size mismatch')
6293 stop
'AVPROD: mismatched vector and matrix'
6324 INTEGER(mpi) :: ispc
6325 INTEGER(mpi) :: item1
6326 INTEGER(mpi) :: item2
6327 INTEGER(mpi) :: itemc
6328 INTEGER(mpi) :: jtem
6329 INTEGER(mpi) :: jtemn
6332 INTEGER(mpi),
INTENT(IN) :: itema
6333 INTEGER(mpi),
INTENT(IN) :: itemb
6334 INTEGER(mpl),
INTENT(OUT) :: ij
6335 INTEGER(mpi),
INTENT(OUT) :: lr
6336 INTEGER(mpi),
INTENT(OUT) :: iprc
6347 item1=max(itema,itemb)
6348 item2=min(itema,itemb)
6349 IF(item2 <= 0.OR.item1 >
napgrp)
RETURN
6352 outer:
DO ispc=1,
nspc
6363 IF(ku < kl) cycle outer
6368 IF(item2 >= jtem.AND.item2 < jtemn)
THEN
6374 IF(item2 < jtem)
THEN
6376 ELSE IF(item2 >= jtemn)
THEN
6391 IF(ku < kl) cycle outer
6395 IF(itemc == jtem)
EXIT
6396 IF(itemc < jtem)
THEN
6398 ELSE IF(itemc > jtem)
THEN
6426 INTEGER(mpi),
INTENT(IN) :: itema
6427 INTEGER(mpi),
INTENT(IN) :: itemb
6452 INTEGER(mpi) :: item1
6453 INTEGER(mpi) :: item2
6454 INTEGER(mpi) :: ipg1
6455 INTEGER(mpi) :: ipg2
6457 INTEGER(mpi) :: iprc
6459 INTEGER(mpi),
INTENT(IN) :: itema
6460 INTEGER(mpi),
INTENT(IN) :: itemb
6462 INTEGER(mpl) ::
ijadd
6466 item1=max(itema,itemb)
6467 item2=min(itema,itemb)
6469 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6470 IF(item1 == item2)
THEN
6479 CALL ijpgrp(ipg1,ipg2,ij,lr,iprc)
6501 INTEGER(mpi) :: item1
6502 INTEGER(mpi) :: item2
6503 INTEGER(mpi) :: jtem
6505 INTEGER(mpi),
INTENT(IN) :: itema
6506 INTEGER(mpi),
INTENT(IN) :: itemb
6515 item1=max(itema,itemb)
6516 item2=min(itema,itemb)
6518 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6526 print *,
' IJCSR3 empty list ', item1, item2, ks, ke
6527 CALL peend(23,
'Aborted, bad matrix index')
6528 stop
'ijcsr3: empty list'
6533 IF(item1 == jtem)
EXIT
6534 IF(item1 < jtem)
THEN
6541 print *,
' IJCSR3 not found ', item1, item2, ks, ke
6542 CALL peend(23,
'Aborted, bad matrix index')
6543 stop
'ijcsr3: not found'
6559 INTEGER(mpi) :: item1
6560 INTEGER(mpi) :: item2
6564 INTEGER(mpl) ::
ijadd
6567 INTEGER(mpi),
INTENT(IN) :: itema
6568 INTEGER(mpi),
INTENT(IN) :: itemb
6573 item1=max(itema,itemb)
6574 item2=min(itema,itemb)
6575 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6584 ij=
ijadd(item1,item2)
6587 ELSE IF (ij < 0)
THEN
6619 INTEGER(mpi) :: ichunk
6623 INTEGER(mpi) :: ispc
6631 INTEGER(mpl) :: ijadd
6653 ir=ipg+(ispc-1)*(
napgrp+1)
6703 REAL(mps),
INTENT(IN) :: deltat
6704 INTEGER(mpi),
INTENT(OUT) :: minut
6705 INTEGER(mpi),
INTENT(OUT):: nhour
6706 REAL(mps),
INTENT(OUT):: secnd
6707 INTEGER(mpi) :: nsecd
6710 nsecd=nint(deltat,
mpi)
6712 minut=nsecd/60-60*nhour
6713 secnd=deltat-60*(minut+60*nhour)
6749 INTEGER(mpi),
INTENT(IN) :: item
6753 INTEGER(mpl) :: length
6754 INTEGER(mpl),
PARAMETER :: four = 4
6758 IF(item <= 0)
RETURN
6779 IF(j == 0)
EXIT inner
6801 IF (
lvllog > 1)
WRITE(
lunlog,*)
'INONE: array increased to', &
6822 INTEGER(mpi) :: iprime
6823 INTEGER(mpi) :: nused
6824 LOGICAL :: finalUpdate
6825 INTEGER(mpl) :: oldLength
6826 INTEGER(mpl) :: newLength
6827 INTEGER(mpl),
PARAMETER :: four = 4
6828 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: tempArr
6829 INTEGER(mpl),
DIMENSION(:),
ALLOCATABLE :: tempVec
6833 IF(finalupdate)
THEN
6842 CALL mpalloc(temparr,four,oldlength,
'INONE: temp array')
6844 CALL mpalloc(tempvec,oldlength,
'INONE: temp vector')
6868 IF(j == 0)
EXIT inner
6869 IF(j == i) cycle outer
6873 IF(.NOT.finalupdate)
RETURN
6878 WRITE(
lunlog,*)
'INONE: array reduced to',newlength,
' words'
6902 IF(j == 0)
EXIT inner
6903 IF(j == i) cycle outer
6920 INTEGER(mpi),
INTENT(IN) :: n
6921 INTEGER(mpi) :: nprime
6922 INTEGER(mpi) :: nsqrt
6927 IF(mod(nprime,2) == 0) nprime=nprime+1
6930 nsqrt=int(sqrt(real(nprime,
mps)),
mpi)
6932 IF(i*(nprime/i) == nprime) cycle outer
6954 INTEGER(mpi) :: idum
6956 INTEGER(mpi) :: indab
6957 INTEGER(mpi) :: itgbi
6958 INTEGER(mpi) :: itgbl
6959 INTEGER(mpi) :: ivgbi
6961 INTEGER(mpi) :: jgrp
6962 INTEGER(mpi) :: lgrp
6964 INTEGER(mpi) :: nc31
6966 INTEGER(mpi) :: nwrd
6967 INTEGER(mpi) :: inone
6972 INTEGER(mpl) :: length
6973 INTEGER(mpl) :: rows
6977 WRITE(
lunlog,*)
'LOOP1: starting'
7000 WRITE(
lunlog,*)
'LOOP1: reading data files'
7011 CALL hmpldf(1,
'Number of words/record in binary file')
7012 CALL hmpdef(8,0.0,60.0,
'not_stored data per record')
7017 CALL peend(20,
'Aborted, bad binary records')
7018 stop
'LOOP1: length of binary record exceeds cache size, wrong file type?'
7032 WRITE(*,*)
'Read all binary data files:'
7051 WRITE(
lunlog,*)
'LOOP1: reading data files again'
7063 CALL peend(21,
'Aborted, no labels/parameters defined')
7064 stop
'LOOP1: no labels/parameters defined'
7069 ' is total number NTGB of labels/parameters'
7071 CALL hmpldf(2,
'Number of entries per label')
7115 WRITE(
lunlog,*)
'LOOP1:',
npresg,
' is number of pre-sigmas'
7116 WRITE(*,*)
'LOOP1:',
npresg,
' is number of pre-sigmas'
7117 IF(
npresg == 0)
WRITE(*,*)
'Warning: no pre-sigmas defined'
7139 WRITE(
lunlog,*)
'LOOP1:',
nvgb,
' is number NVGB of variable parameters'
7144 WRITE(
lunlog,*)
'LOOP1: counting records, NO iteration of entries cut !'
7150 CALL hmpdef(15,0.0,120.0,
'Number of parameters per group')
7185 ' is total number NTPGRP of label/parameter groups'
7201 WRITE(*,112)
' Default pre-sigma =',
regpre, &
7202 ' (if no individual pre-sigma defined)'
7203 WRITE(*,*)
'Pre-sigma factor is',
regula
7206 WRITE(*,*)
'No regularization will be done'
7208 WRITE(*,*)
'Regularization will be done, using factor',
regula
7212 CALL peend(22,
'Aborted, no variable global parameters')
7213 stop
'... no variable global parameters'
7220 IF(presg > 0.0)
THEN
7222 ELSE IF(presg == 0.0.AND.
regpre > 0.0)
THEN
7223 prewt=1.0/real(
regpre**2,mpd)
7239 WRITE(*,101)
'NTGB',
ntgb,
'total number of parameters'
7240 WRITE(*,101)
'NVGB',
nvgb,
'number of variable parameters'
7244 WRITE(*,101)
'Too many variable parameters for diagonalization, fallback is inversion'
7252 WRITE(*,101)
' NREC',
nrec,
'number of records'
7253 IF (
nrecd > 0)
WRITE(*,101)
' NRECD',
nrec,
'number of records containing doubles'
7254 WRITE(*,101)
' NEQN',
neqn,
'number of equations (measurements)'
7255 WRITE(*,101)
' NEGB',
negb,
'number of equations with global parameters'
7256 WRITE(*,101)
' NDGB',
ndgb,
'number of global derivatives'
7258 WRITE(*,101)
' NZGB',
nzgb,
'number of zero global der. (ignored in entry counts)'
7261 WRITE(*,101)
'MREQENF',
mreqenf,
'required number of entries (eqns in binary files)'
7263 WRITE(*,101)
'MREQENF',
mreqenf,
'required number of entries (recs in binary files)'
7266 WRITE(*,101)
'ITEREN',
iteren,
'iterate cut for parameters with less entries'
7267 WRITE(*,101)
'MREQENA',
mreqena,
'required number of entries (from accepted fits)'
7268 IF (
mreqpe > 1)
WRITE(*,101) &
7269 'MREQPE',
mreqpe,
'required number of pair entries'
7270 IF (
msngpe >= 1)
WRITE(*,101) &
7271 'MSNGPE',
msngpe,
'max pair entries single prec. storage'
7272 WRITE(*,101)
'NTGB',
ntgb,
'total number of parameters'
7273 WRITE(*,101)
'NVGB',
nvgb,
'number of variable parameters'
7276 WRITE(*,*)
'Global parameter labels:'
7283 mqi=((mqi-20)/20)*20+1
7291 WRITE(8,101)
' NREC',
nrec,
'number of records'
7292 IF (
nrecd > 0)
WRITE(8,101)
' NRECD',
nrec,
'number of records containing doubles'
7293 WRITE(8,101)
' NEQN',
neqn,
'number of equations (measurements)'
7294 WRITE(8,101)
' NEGB',
negb,
'number of equations with global parameters'
7295 WRITE(8,101)
' NDGB',
ndgb,
'number of global derivatives'
7297 WRITE(8,101)
'MREQENF',
mreqenf,
'required number of entries (eqns in binary files)'
7299 WRITE(8,101)
'MREQENF',
mreqenf,
'required number of entries (recs in binary files)'
7302 WRITE(8,101)
'ITEREN',
iteren,
'iterate cut for parameters with less entries'
7303 WRITE(8,101)
'MREQENA',
mreqena,
'required number of entries (from accepted fits)'
7305 WRITE(
lunlog,*)
'LOOP1: ending'
7309101
FORMAT(1x,a8,
' =',i14,
' = ',a)
7325 INTEGER(mpi) :: ibuf
7327 INTEGER(mpi) :: indab
7333 INTEGER(mpi) :: nc31
7335 INTEGER(mpi) :: nlow
7337 INTEGER(mpi) :: nwrd
7339 INTEGER(mpl) :: length
7340 INTEGER(mpl),
DIMENSION(:),
ALLOCATABLE :: newCounter
7345 WRITE(
lunlog,*)
'LOOP1: iterating'
7347 WRITE(*,*)
'LOOP1: iterating'
7350 CALL mpalloc(newcounter,length,
'new entries counter')
7374 CALL isjajb(nst,ist,ja,jb,jsp)
7375 IF(ja == 0.AND.jb == 0)
EXIT
7381 IF(ij == -2) nlow=nlow+1
7386 newcounter(ij)=newcounter(ij)+1
7415 WRITE(
lunlog,*)
'LOOP1:',
nvgb,
' is number NVGB of variable parameters'
7445 INTEGER(mpi) :: ibuf
7446 INTEGER(mpi) :: icblst
7447 INTEGER(mpi) :: icboff
7448 INTEGER(mpi) :: icgb
7449 INTEGER(mpi) :: icgrp
7450 INTEGER(mpi) :: icount
7451 INTEGER(mpi) :: iext
7452 INTEGER(mpi) :: ihis
7456 INTEGER(mpi) :: ioff
7457 INTEGER(mpi) :: ipoff
7458 INTEGER(mpi) :: iproc
7459 INTEGER(mpi) :: irecmm
7461 INTEGER(mpi) :: itgbi
7462 INTEGER(mpi) :: itgbij
7463 INTEGER(mpi) :: itgbik
7464 INTEGER(mpi) :: ivgbij
7465 INTEGER(mpi) :: ivgbik
7466 INTEGER(mpi) :: ivpgrp
7470 INTEGER(mpi) :: jcgrp
7471 INTEGER(mpi) :: jext
7472 INTEGER(mpi) :: jcgb
7473 INTEGER(mpi) :: jrec
7475 INTEGER(mpi) :: joff
7477 INTEGER(mpi) :: kcgrp
7478 INTEGER(mpi) :: kfile
7480 INTEGER(mpi) :: label
7481 INTEGER(mpi) :: labelf
7482 INTEGER(mpi) :: labell
7483 INTEGER(mpi) :: lvpgrp
7486 INTEGER(mpi) :: maeqnf
7487 INTEGER(mpi) :: nall
7488 INTEGER(mpi) :: naeqna
7489 INTEGER(mpi) :: naeqnf
7490 INTEGER(mpi) :: naeqng
7491 INTEGER(mpi) :: npdblk
7492 INTEGER(mpi) :: nc31
7493 INTEGER(mpi) :: ncachd
7494 INTEGER(mpi) :: ncachi
7495 INTEGER(mpi) :: ncachr
7496 INTEGER(mpi) :: ncon
7499 INTEGER(mpi) :: ndfmax
7500 INTEGER(mpi) :: nfixed
7501 INTEGER(mpi) :: nggd
7502 INTEGER(mpi) :: nggi
7503 INTEGER(mpi) :: nmatmo
7504 INTEGER(mpi) :: noff
7505 INTEGER(mpi) :: npair
7506 INTEGER(mpi) :: npar
7507 INTEGER(mpi) :: nparmx
7509 INTEGER(mpi) :: nrece
7510 INTEGER(mpi) :: nrecf
7511 INTEGER(mpi) :: nrecmm
7513 INTEGER(mpi) :: nwrd
7514 INTEGER(mpi) :: inone
7523 INTEGER(mpl):: nblock
7524 INTEGER(mpl):: nbwrds
7525 INTEGER(mpl):: noff8
7526 INTEGER(mpl):: ndimbi
7527 INTEGER(mpl):: ndimsa(4)
7529 INTEGER(mpl):: nnzero
7530 INTEGER(mpl):: matsiz(2)
7531 INTEGER(mpl):: matwords
7532 INTEGER(mpl):: mbwrds
7533 INTEGER(mpl):: length
7536 INTEGER(mpl),
PARAMETER :: two=2
7537 INTEGER(mpi) :: maxGlobalPar = 0
7538 INTEGER(mpi) :: maxLocalPar = 0
7539 INTEGER(mpi) :: maxEquations = 0
7541 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsGroupList
7542 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsGroupIndex
7543 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecPairedParGroups
7544 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecBlockCounts
7547 SUBROUTINE ndbits(npgrp,ndims,nsparr,ihst)
7549 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7550 INTEGER(mpl),
DIMENSION(4),
INTENT(OUT) :: ndims
7551 INTEGER(mpl),
DIMENSION(:,:),
INTENT(OUT) :: nsparr
7552 INTEGER(mpi),
INTENT(IN) :: ihst
7554 SUBROUTINE ckbits(npgrp,ndims)
7556 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7557 INTEGER(mpl),
DIMENSION(4),
INTENT(OUT) :: ndims
7559 SUBROUTINE spbits(npgrp,nsparr,nsparc)
7561 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7562 INTEGER(mpl),
DIMENSION(:,:),
INTENT(IN) :: nsparr
7563 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: nsparc
7565 SUBROUTINE gpbmap(ngroup,npgrp,npair)
7567 INTEGER(mpi),
INTENT(IN) :: ngroup
7568 INTEGER(mpi),
DIMENSION(:,:),
INTENT(IN) :: npgrp
7569 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npair
7571 SUBROUTINE ggbmap(ipgrp,npair,npgrp)
7573 INTEGER(mpi),
INTENT(IN) :: ipgrp
7574 INTEGER(mpi),
INTENT(OUT) :: npair
7575 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npgrp
7577 SUBROUTINE pbsbits(npgrp,ibsize,nnzero,nblock,nbkrow)
7579 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7580 INTEGER(mpi),
INTENT(IN) :: ibsize
7581 INTEGER(mpl),
INTENT(OUT) :: nnzero
7582 INTEGER(mpl),
INTENT(OUT) :: nblock
7583 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: nbkrow
7585 SUBROUTINE pblbits(npgrp,ibsize,nsparr,nsparc)
7587 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7588 INTEGER(mpi),
INTENT(IN) :: ibsize
7589 INTEGER(mpl),
DIMENSION(:),
INTENT(IN) :: nsparr
7590 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparc
7592 SUBROUTINE prbits(npgrp,nsparr)
7594 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7595 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparr
7597 SUBROUTINE pcbits(npgrp,nsparr,nsparc)
7599 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7600 INTEGER(mpl),
DIMENSION(:),
INTENT(IN) :: nsparr
7601 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparc
7611 WRITE(
lunlog,*)
'LOOP2: starting'
7634 WRITE(
lunlog,*)
' Elimination for constraints enforced for solution by decomposition!'
7639 WRITE(
lunlog,*)
' Lagrange multipliers enforced for solution by sparsePARDISO!'
7644 WRITE(
lunlog,*)
' Elimination for constraints with mpqldec enforced (LAPACK only for unpacked storage)!'
7661 noff8=int(
nagb,mpl)*int(
nagb-1,mpl)/2
7697 print *,
' Variable parameter groups ',
nvpgrp
7745 CALL mpalloc(vecconsgrouplist,length,
'constraint group list')
7746 CALL mpalloc(vecconsgroupindex,length,
'constraint group index')
7756 WRITE(
lunlog,*)
'LOOP2: start event reading'
7796 WRITE(*,*)
'Record number ',
nrec,
' from file ',kfile
7797 IF (wgh /= 1.0)
WRITE(*,*)
' weight ',wrec
7801 CALL isjajb(nst,ist,ja,jb,jsp)
7805 IF(nda ==
mdebg2+1)
WRITE(*,*)
'... and more data'
7810 WRITE(*,*)
'Local derivatives:'
7812107
FORMAT(6(i3,g12.4))
7814 WRITE(*,*)
'Global derivatives:'
7817108
FORMAT(3i11,g12.4)
7820 WRITE(*,*)
'total_par_label __label__ var_par_index derivative'
7835 CALL isjajb(nst,ist,ja,jb,jsp)
7836 IF(ja == 0.AND.jb == 0)
EXIT
7876 IF (kcgrp == jcgrp) cycle
7887 IF (vecconsgroupindex(k) == 0)
THEN
7890 vecconsgrouplist(icgrp)=k
7905 IF (vecconsgroupindex(k) < icount)
THEN
7907 vecconsgroupindex(k)=icount
7925 IF (nfixed > 0) naeqnf=naeqnf+1
7928 IF(ja /= 0.AND.jb /= 0)
THEN
7937 IF (naeqnf > maeqnf) nrecf=nrecf+1
7941 maxglobalpar=max(
nagbn,maxglobalpar)
7942 maxlocalpar=max(
nalcn,maxlocalpar)
7943 maxequations=max(
naeqn,maxequations)
7946 dstat(1)=dstat(1)+real((nwrd+2)*2,mpd)
7947 dstat(2)=dstat(2)+real(
nagbn+2,mpd)
7952 vecconsgroupindex(vecconsgrouplist(k))=0
7957 IF (
nagbn == 0)
THEN
7976 IF (ivpgrp /= lvpgrp)
THEN
8031 (irecmm >= nrecmm.OR.irecmm ==
mxrec))
THEN
8032 IF (nmatmo == 0)
THEN
8034 WRITE(*,*)
'Monitoring of sparse matrix construction'
8035 WRITE(*,*)
' records ........ off-diagonal elements ', &
8036 '....... compression memory'
8037 WRITE(*,*)
' non-zero used(double) used', &
8042 gbc=1.0e-9*real((mpi*ndimsa(2)+mpd*ndimsa(3)+mps*ndimsa(4))/mpi*(bit_size(1_mpi)/8),mps)
8043 gbu=1.0e-9*real(((mpi+mpd)*(ndimsa(3)+ndimsa(4)))/mpi*(bit_size(1_mpi)/8),mps)
8045 WRITE(*,1177) irecmm,ndimsa(1),ndimsa(3),ndimsa(4),cpr,gbc
80461177
FORMAT(i9,3i13,f10.2,f11.6)
8047 DO WHILE(irecmm >= nrecmm)
8066 WRITE(
lunlog,*)
'LOOP2: event reading ended - end of data'
8068 dstat(k)=dstat(k)/real(
nrec,mpd)
8082 CALL mpalloc(vecpairedpargroups,length,
'paired global parameter groups (I)')
8084 print *,
' Total parameter groups pairs',
ntpgrp
8087 CALL ggbmap(i,npair,vecpairedpargroups)
8155 IF(ivgbij > 0.AND.ivgbik > 0)
THEN
8157 IF (
mprint > 1)
WRITE(*,*)
'add index pair ',ivgbij,ivgbik
8234 nparmx=max(nparmx,int(rows,mpi))
8256 WRITE(*,*)
' Parameter block', i, ib-ia, jb-ja, labelf, labell
8260 WRITE(
lunlog,*)
'Detected',
npblck,
'(disjoint) parameter blocks, max size ', nparmx
8262 WRITE(*,*)
'Detected',
npblck,
'(disjoint) parameter blocks, max size ', nparmx
8264 WRITE(*,*)
'Using block diagonal storage mode'
8279 IF (
nagb >= 65536)
THEN
8280 noff=int(noff8/1000,mpi)
8290 CALL hmpdef(ihis,0.0,real(
mhispe,mps),
'NDBITS: #off-diagonal elements')
8295 ndgn=ndimsa(3)+ndimsa(4)
8296 matwords=ndimsa(2)+length*4
8311 length=int(npdblk,mpl)
8312 CALL mpalloc(vecblockcounts,length,
'sparse matrix row offsets (CSR3)')
8314 nbwrds=2*int(nblock,mpl)*int(
ipdbsz(i)*
ipdbsz(i)+1,mpl)
8315 IF ((i == 1).OR.(nbwrds < mbwrds))
THEN
8339 IF (
fcache(2) == 0.0)
THEN
8341 fcache(2)=real(dstat(2),mps)
8342 fcache(3)=real(dstat(3),mps)
8363 ncachd=
ncache-ncachr-ncachi
8394 WRITE(lu,101)
'NTGB',
ntgb,
'total number of parameters'
8395 WRITE(lu,102)
'(all parameters, appearing in binary files)'
8396 WRITE(lu,101)
'NVGB',
nvgb,
'number of variable parameters'
8397 WRITE(lu,102)
'(appearing in fit matrix/vectors)'
8398 WRITE(lu,101)
'NAGB',
nagb,
'number of all parameters'
8399 WRITE(lu,102)
'(including Lagrange multiplier or reduced)'
8400 WRITE(lu,101)
'NTPGRP',
ntpgrp,
'total number of parameter groups'
8401 WRITE(lu,101)
'NVPGRP',
nvpgrp,
'number of variable parameter groups'
8402 WRITE(lu,101)
'NFGB',
nfgb,
'number of fit parameters'
8404 WRITE(lu,101)
'MBANDW',
mbandw,
'band width of preconditioner matrix'
8405 WRITE(lu,102)
'(if <0, no preconditioner matrix)'
8407 IF (
nagb >= 65536)
THEN
8408 WRITE(lu,101)
'NOFF/K',noff,
'max number of off-diagonal elements'
8410 WRITE(lu,101)
'NOFF',noff,
'max number of off-diagonal elements'
8413 IF (
nagb >= 65536)
THEN
8414 WRITE(lu,101)
'NDGN/K',ndgn/1000,
'actual number of off-diagonal elements'
8416 WRITE(lu,101)
'NDGN',ndgn,
'actual number of off-diagonal elements'
8419 WRITE(lu,101)
'NCGB',
ncgb,
'number of constraints'
8420 WRITE(lu,101)
'NAGBN',
nagbn,
'max number of global parameters in an event'
8421 WRITE(lu,101)
'NALCN',
nalcn,
'max number of local parameters in an event'
8422 WRITE(lu,101)
'NAEQN',
naeqn,
'max number of equations in an event'
8424 WRITE(lu,101)
'NAEQNA',naeqna,
'number of equations'
8425 WRITE(lu,101)
'NAEQNG',naeqng, &
8426 'number of equations with global parameters'
8427 WRITE(lu,101)
'NAEQNF',naeqnf, &
8428 'number of equations with fixed global parameters'
8429 WRITE(lu,101)
'NRECF',nrecf, &
8430 'number of records with fixed global parameters'
8433 WRITE(lu,101)
'NRECE',nrece, &
8434 'number of records without variable parameters'
8437 WRITE(lu,101)
'NCACHE',
ncache,
'number of words for caching'
8438 WRITE(lu,111) (
fcache(k)*100.0,k=1,3)
8439111
FORMAT(22x,
'cache splitting ',3(f6.1,
' %'))
8444 WRITE(lu,*)
'Solution method and matrix-storage mode:'
8446 WRITE(lu,*)
' METSOL = 1: matrix inversion'
8447 ELSE IF(
metsol == 2)
THEN
8448 WRITE(lu,*)
' METSOL = 2: diagonalization'
8449 ELSE IF(
metsol == 3)
THEN
8450 WRITE(lu,*)
' METSOL = 3: decomposition'
8451 ELSE IF(
metsol == 4)
THEN
8452 WRITE(lu,*)
' METSOL = 4: MINRES (rtol',
mrestl,
')'
8453 ELSE IF(
metsol == 5)
THEN
8454 WRITE(lu,*)
' METSOL = 5: MINRES-QLP (rtol',
mrestl,
')'
8455 ELSE IF(
metsol == 6)
THEN
8456 WRITE(lu,*)
' METSOL = 6: GMRES'
8458 ELSE IF(
metsol == 7)
THEN
8459 WRITE(lu,*)
' METSOL = 7: LAPACK factorization'
8460 ELSE IF(
metsol == 8)
THEN
8461 WRITE(lu,*)
' METSOL = 8: LAPACK factorization'
8463 ELSE IF(
metsol == 9)
THEN
8464 WRITE(lu,*)
' METSOL = 9: Intel oneMKL PARDISO'
8468 WRITE(lu,*)
' with',
mitera,
' iterations'
8470 WRITE(lu,*)
' MATSTO = 0: unpacked symmetric matrix, ',
'n*n elements'
8471 ELSE IF(
matsto == 1)
THEN
8472 WRITE(lu,*)
' MATSTO = 1: full symmetric matrix, ',
'(n*n+n)/2 elements'
8473 ELSE IF(
matsto == 2)
THEN
8474 WRITE(lu,*)
' MATSTO = 2: sparse matrix (custom)'
8475 ELSE IF(
matsto == 3)
THEN
8477 WRITE(lu,*)
' MATSTO = 3: sparse matrix (upper triangle, CSR3)'
8479 WRITE(lu,*)
' MATSTO = 3: sparse matrix (upper triangle, BSR3)'
8480 WRITE(lu,*)
' block size',
matbsz
8484 WRITE(lu,*)
' block diagonal with',
npblck,
' blocks'
8486 IF(
mextnd>0)
WRITE(lu,*)
' with extended storage'
8487 IF(
dflim /= 0.0)
THEN
8488 WRITE(lu,103)
'Convergence assumed, if expected dF <',
dflim
8493 WRITE(lu,*)
'Constraints handled by elimination with LAPACK'
8495 WRITE(lu,*)
'Constraints handled by elimination'
8498 WRITE(lu,*)
'Constraints handled by Lagrange multipliers'
8505 CALL peend(28,
'Aborted, no local parameters')
8506 stop
'LOOP2: stopping due to missing local parameters'
8527105
FORMAT(
' Constants C1, C2 for Wolfe conditions:',g12.4,
', ',g12.4)
8535 matwords=(length+
nagb+1)*2
8542 ELSE IF (
matsto == 2)
THEN
8543 matsiz(1)=ndimsa(3)+
nagb
8558 IF (icblst > icboff)
THEN
8564 nall = npar;
IF (
icelim <= 0) nall=npar+ncon
8568 matsiz(1)=matsiz(1)+k
8570 matsiz(1)=matsiz(1)+nall
8575 matwords=matwords+matsiz(1)*2+matsiz(2)
8581 IF (matwords < 250000)
THEN
8582 WRITE(*,*)
'Size of global matrix: < 1 MB'
8584 WRITE(*,*)
'Size of global matrix:',int(real(matwords,mps)*4.0e-6,mpi),
' MB'
8590 WRITE(
lunlog,*)
' Cut values of Chi^2/Ndf and Chi2,'
8591 WRITE(
lunlog,*)
' corresponding to 2 and 3 standard deviations'
8592 WRITE(
lunlog,*)
' Ndf Chi^2/Ndf(2) Chi^2(2) ', &
8593 ' Chi^2/Ndf(3) Chi^2(3)'
8596 IF(ndf >
naeqn)
EXIT
8599 ELSE IF(ndf < 20)
THEN
8601 ELSE IF(ndf < 100)
THEN
8603 ELSE IF(ndf < 200)
THEN
8610 WRITE(
lunlog,106) ndf,chin2,chin2*real(ndf,mps),chin3, chin3*real(ndf,mps)
8613 WRITE(
lunlog,*)
'LOOP2: ending'
8617 IF (
ncgbe /= 0)
THEN
8620 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
8621 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
8622 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
8623 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
8624 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
8625 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
8626 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
8628 WRITE(*,*)
' Number of empty constraints =',abs(
ncgbe),
', should be 0'
8629 WRITE(*,*)
' => please check constraint definition, mille data'
8631 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
8632 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
8633 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
8634 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
8635 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
8636 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
8637 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
8642101
FORMAT(1x,a8,
' =',i14,
' = ',a)
8644103
FORMAT(1x,a,g12.4)
8645106
FORMAT(i6,2(3x,f9.3,f12.1,3x))
8660 INTEGER(mpi) :: imed
8663 INTEGER(mpi) :: nent
8664 INTEGER(mpi),
DIMENSION(measBins) :: isuml
8665 INTEGER(mpi),
DIMENSION(measBins) :: isums
8669 INTEGER(mpl) :: ioff
8672 DATA lfirst /.true./
8685 WRITE(
lunmon,
'(A)')
'*** Normalized residuals grouped by first global label (per local fit cycle) ***'
8687 WRITE(
lunmon,
'(A)')
'*** Pulls grouped by first global label (per local fit cycle) ***'
8689 WRITE(
lunmon,
'(A)')
'! LFC Label Entries Median RMS(MAD) <error>'
8694#ifdef SCOREP_USER_ENABLE
8695 scorep_user_region_by_name_begin(
"UR_monres", scorep_user_region_type_common)
8711 IF (2*isuml(j) > nent)
EXIT
8715 IF (isuml(j) > isuml(j-1)) amed=amed+real(nent-2*isuml(j-1),mps)/real(2*isuml(j)-2*isuml(j-1),mps)
8728 isums(j)=isums(j)+isums(j-1)
8732 IF (2*isums(j) > nent)
EXIT
8735 IF (isums(j) > isums(j-1)) amad=amad+real(nent-2*isums(j-1),mps)/real(2*isums(j)-2*isums(j-1),mps)
8747#ifdef SCOREP_USER_ENABLE
8748 scorep_user_region_by_name_end(
"UR_monres")
8752110
FORMAT(i5,2i10,3g14.5)
8767 INTEGER(mpi) :: ioff
8768 INTEGER(mpi) :: ipar0
8769 INTEGER(mpi) :: ncon
8770 INTEGER(mpi) :: npar
8771 INTEGER(mpi) :: nextra
8773 INTEGER :: nbopt, nboptx, ILAENV
8776 INTEGER(mpl),
INTENT(IN) :: msize(2)
8778 INTEGER(mpl) :: length
8779 INTEGER(mpl) :: nwrdpc
8780 INTEGER(mpl),
PARAMETER :: three = 3
8793 CALL mpalloc(
aux,length,
' local fit scratch array: aux')
8794 CALL mpalloc(
vbnd,length,
' local fit scratch array: vbnd')
8795 CALL mpalloc(
vbdr,length,
' local fit scratch array: vbdr')
8798 CALL mpalloc(
vbk,length,
' local fit scratch array: vbk')
8801 CALL mpalloc(
vzru,length,
' local fit scratch array: vzru')
8802 CALL mpalloc(
scdiag,length,
' local fit scratch array: scdiag')
8803 CALL mpalloc(
scflag,length,
' local fit scratch array: scflag')
8804 CALL mpalloc(
ibandh,2*length,
' local fit band width hist.: ibandh')
8818 length=4*
ncblck;
IF(ncon == 0) length=0
8827 nwrdpc=nwrdpc+length
8839 length=(ncon*ncon+ncon)/2
8864 length=length+npar+nextra
8870 IF(
mbandw == 0) length=length+1
8878 nwrdpc=nwrdpc+2*length
8879 IF (nwrdpc > 250000)
THEN
8881 WRITE(*,*)
'Size of preconditioner matrix:',int(real(nwrdpc,mps)*4.0e-6,mpi),
' MB'
8903 CALL peend(23,
'Aborted, bad matrix index (will exceed 32bit)')
8904 stop
'vmprep: bad index (matrix to large for diagonalization)'
8926 nbopt = ilaenv( 1_mpl,
'DSYTRF',
'U', int(
nagb,mpl), int(
nagb,mpl), -1_mpl, -1_mpl )
8928 print *,
'LAPACK optimal block size for DSYTRF:', nbopt
8929 lplwrk=length*int(nbopt,mpl)
8937 nbopt = ilaenv( 1_mpl,
'DORMQL',
'RN', int(npar,mpl), int(npar,mpl), int(ncon,mpl), int(npar,mpl) )
8938 IF (int(npar,mpl)*int(nbopt,mpl) >
lplwrk)
THEN
8939 lplwrk=int(npar,mpl)*int(nbopt,mpl)
8944 print *,
'LAPACK optimal block size for DORMQL:', nboptx
8963 INTEGER(mpi) :: icoff
8964 INTEGER(mpi) :: ipoff
8967 INTEGER(mpi) :: ncon
8968 INTEGER(mpi) :: nfit
8969 INTEGER(mpi) :: npar
8970 INTEGER(mpi) :: nrank
8971 INTEGER(mpl) :: imoff
8972 INTEGER(mpl) :: ioff1
8990 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9005 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9007 IF(nfit < npar)
THEN
9025 WRITE(
lunlog,*)
'Inversion of global matrix (A->A^-1)'
9032 IF(nfit /= nrank)
THEN
9033 WRITE(*,*)
'Warning: the rank defect of the symmetric',nfit, &
9034 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9035 WRITE(lun,*)
'Warning: the rank defect of the symmetric',nfit, &
9036 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9039 WRITE(*,*)
' --> enforcing SUBITO mode'
9040 WRITE(lun,*)
' --> enforcing SUBITO mode'
9042 ELSE IF(
ndefec == 0)
THEN
9044 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9046 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9057 IF(nfit < npar)
THEN
9076 INTEGER(mpi) :: icoff
9077 INTEGER(mpi) :: ipoff
9080 INTEGER(mpi) :: ncon
9081 INTEGER(mpi) :: nfit
9082 INTEGER(mpi) :: npar
9083 INTEGER(mpi) :: nrank
9084 INTEGER(mpl) :: imoff
9085 INTEGER(mpl) :: ioff1
9100 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9114 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9116 IF(nfit < npar)
THEN
9134 WRITE(
lunlog,*)
'Decomposition of global matrix (A->L*D*L^t)'
9140 IF(nfit /= nrank)
THEN
9141 WRITE(*,*)
'Warning: the rank defect of the symmetric',nfit, &
9142 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9143 WRITE(lun,*)
'Warning: the rank defect of the symmetric',nfit, &
9144 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9147 WRITE(*,*)
' --> enforcing SUBITO mode'
9148 WRITE(lun,*)
' --> enforcing SUBITO mode'
9150 ELSE IF(
ndefec == 0)
THEN
9152 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9154 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9156 WRITE(lun,*)
' largest diagonal element (LDLt)', evmax
9157 WRITE(lun,*)
' smallest diagonal element (LDLt)', evmin
9166 IF(nfit < npar)
THEN
9188 INTEGER(mpi) :: icoff
9189 INTEGER(mpi) :: ipoff
9192 INTEGER(mpi) :: ncon
9193 INTEGER(mpi) :: nfit
9194 INTEGER(mpi) :: npar
9195 INTEGER(mpl) :: imoff
9196 INTEGER(mpl) :: ioff1
9197 INTEGER(mpi) :: infolp
9217 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9232 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9234 IF(nfit < npar)
THEN
9251 IF (nfit > npar)
THEN
9254 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*D*L^t)'
9258#ifdef SCOREP_USER_ENABLE
9259 scorep_user_region_by_name_begin(
"UR_dsptrf", scorep_user_region_type_common)
9262#ifdef SCOREP_USER_ENABLE
9263 scorep_user_region_by_name_end(
"UR_dsptrf")
9270 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*L^t)'
9274#ifdef SCOREP_USER_ENABLE
9275 scorep_user_region_by_name_begin(
"UR_dpptrf", scorep_user_region_type_common)
9277 CALL dpptrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),infolp)
9278#ifdef SCOREP_USER_ENABLE
9279 scorep_user_region_by_name_end(
"UR_dpptrf")
9287 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9289 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9293 WRITE(*,*)
'Warning: factorization of the symmetric',nfit, &
9294 '-by-',nfit,
' failed at index ', infolp
9295 WRITE(lun,*)
'Warning: factorization of the symmetric',nfit, &
9296 '-by-',nfit,
' failed at index ', infolp
9297 CALL peend(29,
'Aborted, factorization of global matrix failed')
9298 stop
'mdptrf: bad matrix'
9303 IF (nfit > npar)
THEN
9306 IF(infolp /= 0) print *,
' DSPTRS failed: ', infolp
9308 CALL dpptrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),&
9310 IF(infolp /= 0) print *,
' DPPTRS failed: ', infolp
9314 IF(nfit < npar)
THEN
9335 INTEGER(mpi) :: icoff
9336 INTEGER(mpi) :: ipoff
9339 INTEGER(mpi) :: ncon
9340 INTEGER(mpi) :: nfit
9341 INTEGER(mpi) :: npar
9342 INTEGER(mpl) :: imoff
9343 INTEGER(mpl) :: ioff1
9344 INTEGER(mpl) :: iloff
9345 INTEGER(mpi) :: infolp
9366 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9386 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9388 IF(nfit < npar)
THEN
9392 CALL dtrtrs(
'L',
'T',
'N',int(ncon,mpl),1_mpl,
lapackql(iloff+npar-ncon+1:),int(npar,mpl),&
9394 IF(infolp /= 0) print *,
' DTRTRS failed: ', infolp
9396 CALL dormql(
'L',
'T',int(npar,mpl),1_mpl,int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9398 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9417 IF (nfit > npar)
THEN
9420 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*D*L^t)'
9424#ifdef SCOREP_USER_ENABLE
9425 scorep_user_region_by_name_begin(
"UR_dsytrf", scorep_user_region_type_common)
9427 CALL dsytrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(nfit,mpl),&
9429#ifdef SCOREP_USER_ENABLE
9430 scorep_user_region_by_name_end(
"UR_dsytrf")
9437 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*L^t)'
9441#ifdef SCOREP_USER_ENABLE
9442 scorep_user_region_by_name_begin(
"UR_dpotrf", scorep_user_region_type_common)
9444 CALL dpotrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(npar,mpl),infolp)
9445#ifdef SCOREP_USER_ENABLE
9446 scorep_user_region_by_name_end(
"UR_dpotrf")
9454 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9456 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9460 WRITE(*,*)
'Warning: factorization of the symmetric',nfit, &
9461 '-by-',nfit,
' failed at index ', infolp
9462 WRITE(lun,*)
'Warning: factorization of the symmetric',nfit, &
9463 '-by-',nfit,
' failed at index ', infolp
9464 CALL peend(29,
'Aborted, factorization of global matrix failed')
9465 stop
'mdutrf: bad matrix'
9470 IF (nfit > npar)
THEN
9471 CALL dsytrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),int(nfit,mpl),&
9473 IF(infolp /= 0) print *,
' DSYTRS failed: ', infolp
9475 CALL dpotrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),int(npar,mpl),&
9477 IF(infolp /= 0) print *,
' DPOTRS failed: ', infolp
9481 IF(nfit < npar)
THEN
9486 CALL dormql(
'L',
'N',int(npar,mpl),1_mpl,int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9488 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9495 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9517 INTEGER(mpi) :: icboff
9518 INTEGER(mpi) :: icblst
9519 INTEGER(mpi) :: icoff
9520 INTEGER(mpi) :: icfrst
9521 INTEGER(mpi) :: iclast
9522 INTEGER(mpi) :: ipfrst
9523 INTEGER(mpi) :: iplast
9524 INTEGER(mpi) :: ipoff
9527 INTEGER(mpi) :: ncon
9528 INTEGER(mpi) :: npar
9530 INTEGER(mpl) :: imoff
9531 INTEGER(mpl) :: iloff
9532 INTEGER(mpi) :: infolp
9533 INTEGER :: nbopt, ILAENV
9535 REAL(mpd),
INTENT(IN) :: a(mszcon)
9536 REAL(mpd),
INTENT(OUT) :: emin
9537 REAL(mpd),
INTENT(OUT) :: emax
9548 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9567 DO icb=icboff+1,icboff+icblst
9574 lapackql(iloff+ipfrst:iloff+iplast)=a(imoff+1:imoff+npb)
9591 nbopt = ilaenv( 1_mpl,
'DGEQLF',
'', int(npar,mpl), int(ncon,mpl), int(npar,mpl), -1_mpl )
9592 print *,
'LAPACK optimal block size for DGEQLF:', nbopt
9593 lplwrk=int(ncon,mpl)*int(nbopt,mpl)
9596#ifdef SCOREP_USER_ENABLE
9597 scorep_user_region_by_name_begin(
"UR_dgeqlf", scorep_user_region_type_common)
9599 CALL dgeqlf(int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9601 IF(infolp /= 0) print *,
' DGEQLF failed: ', infolp
9602#ifdef SCOREP_USER_ENABLE
9603 scorep_user_region_by_name_end(
"UR_dgeqlf")
9607 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9610 IF(emax < emin)
THEN
9638 INTEGER(mpi) :: icoff
9639 INTEGER(mpi) :: ipoff
9641 INTEGER(mpi) :: ncon
9642 INTEGER(mpi) :: npar
9643 INTEGER(mpl) :: imoff
9644 INTEGER(mpl) :: iloff
9645 INTEGER(mpi) :: infolp
9646 CHARACTER (LEN=1) :: transr, transl
9648 LOGICAL,
INTENT(IN) :: t
9667 IF(ncon <= 0 ) cycle
9670#ifdef SCOREP_USER_ENABLE
9671 scorep_user_region_by_name_begin(
"UR_dormql", scorep_user_region_type_common)
9679 DO i=ipoff+1,ipoff+npar
9685 CALL dormql(
'R',transr,int(npar,mpl),int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),&
9688 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9690 CALL dormql(
'L',transl,int(npar,mpl),int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),&
9693 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9694#ifdef SCOREP_USER_ENABLE
9695 scorep_user_region_by_name_end(
"UR_dormql")
9699 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9705include
'mkl_pardiso.f90'
9736 TYPE(mkl_pardiso_handle) :: pt(64)
9738 INTEGER(mpl),
PARAMETER :: maxfct =1
9739 INTEGER(mpl),
PARAMETER :: mnum = 1
9740 INTEGER(mpl),
PARAMETER :: nrhs = 1
9742 INTEGER(mpl) :: mtype
9743 INTEGER(mpl) :: phase
9744 INTEGER(mpl) :: error
9745 INTEGER(mpl) :: msglvl
9749 INTEGER(mpl) :: idum(1)
9751 INTEGER(mpl) :: length
9752 INTEGER(mpi) :: nfill
9753 INTEGER(mpi) :: npdblk
9754 REAL(mpd) :: adum(1)
9755 REAL(mpd) :: ddum(1)
9757 INTEGER(mpl) :: iparm(64)
9758 REAL(mpd),
ALLOCATABLE :: b( : )
9759 REAL(mpd),
ALLOCATABLE :: x( : )
9773#ifdef SCOREP_USER_ENABLE
9774 scorep_user_region_by_name_begin(
"UR_mspd00", scorep_user_region_type_common)
9781 WRITE(*,*)
'MSPARDISO: number of rows to fill up = ', nfill
9794 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl), adum, idum, idum, &
9795 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9796 IF (error /= 0)
THEN
9797 WRITE(lun,*)
'The following ERROR was detected: ', error
9798 WRITE(*,
'(A,2I10)')
' PARDISO release failed (phase, error): ', phase, error
9799 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9800 CALL peend(40,
'Aborted, other error: PARDISO release')
9801 stop
'MSPARDISO: stopping due to error in PARDISO release'
9818 IF (iparm(1) == 0)
WRITE(lun,*)
'PARDISO using defaults '
9819 IF (iparm(43) /= 0)
THEN
9820 WRITE(lun,*)
'PARDISO: computation of the diagonal of inverse matrix not implemented !'
9828#ifdef SCOREP_USER_ENABLE
9829 scorep_user_region_by_name_end(
"UR_mspd00")
9837 WRITE(
lunlog,*)
'Decomposition of global matrix (A->L*D*L^t)'
9844#ifdef SCOREP_USER_ENABLE
9845 scorep_user_region_by_name_begin(
"UR_mspd11", scorep_user_region_type_common)
9854 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9857 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9858 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9859#ifdef SCOREP_USER_ENABLE
9860 scorep_user_region_by_name_end(
"UR_mspd11")
9863 WRITE(lun,*)
'PARDISO reordering completed ... '
9864 WRITE(lun,*)
'PARDISO peak memory required (KB)', iparm(15)
9867 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9870 IF (error /= 0)
THEN
9871 WRITE(lun,*)
'The following ERROR was detected: ', error
9872 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9873 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9874 CALL peend(40,
'Aborted, other error: PARDISO reordering')
9875 stop
'MSPARDISO: stopping due to error in PARDISO reordering'
9877 IF (iparm(60) == 0)
THEN
9882 WRITE(lun,*)
'Size (KB) of allocated memory = ',
ipdmem
9883 WRITE(lun,*)
'Number of nonzeros in factors = ',iparm(18)
9884 WRITE(lun,*)
'Number of factorization MFLOPS = ',iparm(19)
9888#ifdef SCOREP_USER_ENABLE
9889 scorep_user_region_by_name_begin(
"UR_mspd22", scorep_user_region_type_common)
9892 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9893 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9894#ifdef SCOREP_USER_ENABLE
9895 scorep_user_region_by_name_end(
"UR_mspd22")
9898 WRITE(lun,*)
'PARDISO factorization completed ... '
9901 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9904 IF (error /= 0)
THEN
9905 WRITE(lun,*)
'The following ERROR was detected: ', error
9906 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9907 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9908 CALL peend(40,
'Aborted, other error: PARDISO factorization')
9909 stop
'MSPARDISO: stopping due to error in PARDISO factorization'
9912 IF (iparm(14) > 0) &
9913 WRITE(lun,*)
'Number of perturbed pivots = ',iparm(14)
9914 WRITE(lun,*)
'Number of positive eigenvalues = ',iparm(22)-nfill
9915 WRITE(lun,*)
'Number of negative eigenvalues = ',iparm(23)
9916 ELSE IF (iparm(30) > 0)
THEN
9917 WRITE(lun,*)
'Equation with bad pivot (<=0.) = ',iparm(30)
9926 CALL mpalloc(b,length,
' PARDISO r.h.s')
9927 CALL mpalloc(x,length,
' PARDISO solution')
9930#ifdef SCOREP_USER_ENABLE
9931 scorep_user_region_by_name_begin(
"UR_mspd33", scorep_user_region_type_common)
9935 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9936 idum, nrhs, iparm, msglvl, b, x, error)
9937#ifdef SCOREP_USER_ENABLE
9938 scorep_user_region_by_name_end(
"UR_mspd33")
9944 WRITE(lun,*)
'PARDISO solve completed ... '
9945 IF (error /= 0)
THEN
9946 WRITE(lun,*)
'The following ERROR was detected: ', error
9947 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9948 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9949 CALL peend(40,
'Aborted, other error: PARDISO solve')
9950 stop
'MSPARDISO: stopping due to error in PARDISO solve'
9964 INTEGER(mpi) :: iast
9965 INTEGER(mpi) :: idia
9966 INTEGER(mpi) :: imin
9967 INTEGER(mpl) :: ioff1
9969 INTEGER(mpi) :: last
9971 INTEGER(mpi) :: nmax
9972 INTEGER(mpi) :: nmin
9973 INTEGER(mpi) :: ntop
9995 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
10032 DO WHILE(ntop < nmax)
10036 CALL hmpdef(7,real(nmin,mps),real(ntop,mps),
'log10 of positive eigenvalues')
10046 iast=max(1,imin-60)
10047 CALL gmpdef(3,2,
'low-value end of eigenvalues')
10050 CALL gmpxy(3,real(i,mps),evalue)
10064 WRITE(lun,*)
'The first (largest) eigenvalues ...'
10067 WRITE(lun,*)
'The last eigenvalues ... up to',last
10071 WRITE(lun,*)
'The eigenvalues from',
nvgb+1,
' to',
nagb
10075 WRITE(lun,*)
'Log10 + 3 of ',
nfgb,
' eigenvalues in decreasing',
' order'
10076 WRITE(lun,*)
'(for Eigenvalue < 0.001 the value 0.0 is shown)'
10079 'printed for negative eigenvalues'
10082 WRITE(lun,*) last,
' significances: insignificant if ', &
10083 'compatible with N(0,1)'
10112 INTEGER(mpl) :: ioff1
10113 INTEGER(mpl) :: ioff2
10149 INTEGER(mpi) :: istop
10150 INTEGER(mpi) :: itn
10151 INTEGER(mpi) :: itnlim
10152 INTEGER(mpi) :: lun
10153 INTEGER(mpi) :: nout
10154 INTEGER(mpi) :: nrkd
10155 INTEGER(mpi) :: nrkd2
10161 REAL(mpd) :: arnorm
10200 WRITE(lun,*)
'MMINRS: PRECONS ended ', nrkd
10204 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10205 ELSE IF(
mbandw > 0)
THEN
10211 WRITE(lun,*)
'MMINRS: EQUDECS ended ', nrkd, nrkd2
10215 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10218 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10232 IF (
istopa == 0) print *,
'MINRES: istop=0, exact solution x=0.'
10247 INTEGER(mpi) :: istop
10248 INTEGER(mpi) :: itn
10249 INTEGER(mpi) :: itnlim
10250 INTEGER(mpi) :: lun
10251 INTEGER(mpi) :: nout
10252 INTEGER(mpi) :: nrkd
10253 INTEGER(mpi) :: nrkd2
10256 REAL(mpd) :: mxxnrm
10257 REAL(mpd) :: trcond
10267 mxxnrm = real(
nagb,mpd)/sqrt(epsilon(mxxnrm))
10269 trcond = 1.0_mpd/epsilon(trcond)
10270 ELSE IF(
mrmode == 2)
THEN
10300 WRITE(lun,*)
'MMINRS: PRECONS ended ', nrkd
10304 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10306 ELSE IF(
mbandw > 0)
THEN
10312 WRITE(lun,*)
'MMINRS: EQUDECS ended ', nrkd, nrkd2
10317 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10321 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10336 IF (
istopa == 3) print *,
'MINRES: istop=0, exact solution x=0.'
10352 INTEGER(mpi),
INTENT(IN) :: n
10353 REAL(mpd),
INTENT(IN) :: x(n)
10354 REAL(mpd),
INTENT(OUT) :: y(n)
10374 INTEGER(mpi),
INTENT(IN) :: n
10375 REAL(mpd),
INTENT(IN) :: x(n)
10376 REAL(mpd),
INTENT(OUT) :: y(n)
10407 REAL(mps) :: concu2
10408 REAL(mps) :: concut
10409 REAL,
DIMENSION(2) :: ta
10412 INTEGER(mpi) :: iact
10413 INTEGER(mpi) :: iagain
10414 INTEGER(mpi) :: idx
10415 INTEGER(mpi) :: info
10417 INTEGER(mpi) :: ipoff
10418 INTEGER(mpi) :: icoff
10419 INTEGER(mpl) :: ioff
10420 INTEGER(mpi) :: itgbi
10421 INTEGER(mpi) :: ivgbi
10422 INTEGER(mpi) :: jcalcm
10424 INTEGER(mpi) :: labelg
10425 INTEGER(mpi) :: litera
10426 INTEGER(mpl) :: lrej
10427 INTEGER(mpi) :: lun
10428 INTEGER(mpi) :: lunp
10429 INTEGER(mpi) :: minf
10430 INTEGER(mpi) :: mrati
10431 INTEGER(mpi) :: nan
10432 INTEGER(mpi) :: ncon
10433 INTEGER(mpi) :: nfaci
10434 INTEGER(mpi) :: nloopsol
10435 INTEGER(mpi) :: npar
10436 INTEGER(mpi) :: nrati
10437 INTEGER(mpl) :: nrej
10438 INTEGER(mpi) :: nsol
10439 INTEGER(mpi) :: inone
10441 INTEGER(mpi) :: infolp
10442 INTEGER(mpi) :: nfit
10443 INTEGER(mpl) :: imoff
10447 REAL(mpd) :: dratio
10448 REAL(mpd) :: dwmean
10457 LOGICAL :: warnerss
10458 LOGICAL :: warners3
10460 CHARACTER (LEN=7) :: cratio
10461 CHARACTER (LEN=7) :: cfacin
10462 CHARACTER (LEN=7) :: crjrat
10473 WRITE(lunp,*)
'Solution algorithm: '
10474 WRITE(lunp,121)
'=================================================== '
10477 WRITE(lunp,121)
'solution method:',
'matrix inversion'
10478 ELSE IF(
metsol == 2)
THEN
10479 WRITE(lunp,121)
'solution method:',
'diagonalization'
10480 ELSE IF(
metsol == 3)
THEN
10481 WRITE(lunp,121)
'solution method:',
'decomposition'
10482 ELSE IF(
metsol == 4)
THEN
10483 WRITE(lunp,121)
'solution method:',
'minres (Paige/Saunders)'
10484 ELSE IF(
metsol == 5)
THEN
10485 WRITE(lunp,121)
'solution method:',
'minres-qlp (Choi/Paige/Saunders)'
10487 WRITE(lunp,121)
' ',
' using QR factorization'
10488 ELSE IF(
mrmode == 2)
THEN
10489 WRITE(lunp,121)
' ',
' using QLP factorization'
10491 WRITE(lunp,121)
' ',
' using QR and QLP factorization'
10492 WRITE(lunp,123)
'transition condition',
mrtcnd
10494 ELSE IF(
metsol == 6)
THEN
10495 WRITE(lunp,121)
'solution method:', &
10496 'gmres (generalized minimzation of residuals)'
10498 ELSE IF(
metsol == 7)
THEN
10500 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DSPTRF)'
10502 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DPPTRF)'
10504 IF(
ilperr == 1)
WRITE(lunp,121)
' ',
'with error calculation (D??TRI)'
10505 ELSE IF(
metsol == 8)
THEN
10507 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DSYTRF)'
10509 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DPOTRF)'
10511 IF(
ilperr == 1)
WRITE(lunp,121)
' ',
'with error calculation (D??TRI)'
10513 ELSE IF(
metsol == 9)
THEN
10515 WRITE(lunp,121)
'solution method:',
'Intel oneMKL PARDISO (sparse matrix (CSR3))'
10517 WRITE(lunp,121)
'solution method:',
'Intel oneMKL PARDISO (sparse matrix (BSR3))'
10522 WRITE(lunp,123)
'convergence limit at Delta F=',
dflim
10523 WRITE(lunp,122)
'maximum number of iterations=',
mitera
10526 WRITE(lunp,122)
'matrix recalculation up to ',
matrit,
'. iteration'
10530 WRITE(lunp,121)
'matrix storage:',
'full'
10531 ELSE IF(
matsto == 2)
THEN
10532 WRITE(lunp,121)
'matrix storage:',
'sparse'
10534 WRITE(lunp,122)
'pre-con band-width parameter=',
mbandw
10536 WRITE(lunp,121)
'pre-conditioning:',
'default'
10537 ELSE IF(
mbandw < 0)
THEN
10538 WRITE(lunp,121)
'pre-conditioning:',
'none!'
10539 ELSE IF(
mbandw > 0)
THEN
10541 WRITE(lunp,121)
'pre-conditioning=',
'skyline-matrix (rank preserving)'
10543 WRITE(lunp,121)
'pre-conditioning=',
'band-matrix'
10548 WRITE(lunp,121)
'using pre-sigmas:',
'no'
10551 WRITE(lunp,124)
'pre-sigmas defined for', &
10552 REAL(100*npresg,mps)/
REAL(nvgb,mps),
' % of variable parameters'
10553 WRITE(lunp,123)
'default pre-sigma=',
regpre
10556 WRITE(lunp,121)
'regularization:',
'no'
10558 WRITE(lunp,121)
'regularization:',
'yes'
10559 WRITE(lunp,123)
'regularization factor=',
regula
10563 WRITE(lunp,121)
'Chi square cut equiv 3 st.dev applied'
10564 WRITE(lunp,123)
'... in first iteration with factor',
chicut
10565 WRITE(lunp,123)
'... in second iteration with factor',
chirem
10566 WRITE(lunp,121)
' (reduced by sqrt in next iterations)'
10569 WRITE(lunp,121)
'Scaling of measurement errors applied'
10570 WRITE(lunp,123)
'... factor for "global" measuements',
dscerr(1)
10571 WRITE(lunp,123)
'... factor for "local" measuements',
dscerr(2)
10574 WRITE(lunp,122)
'Down-weighting of outliers in',
lhuber,
' iterations'
10575 WRITE(lunp,123)
'Cut on downweight fraction',
dwcut
10579121
FORMAT(1x,a40,3x,a)
10580122
FORMAT(1x,a40,3x,i0,a)
10581123
FORMAT(1x,a40,2x,e9.2)
10582124
FORMAT(1x,a40,3x,f5.1,a)
10592 stepl =real(stp,mps)
10604 ELSE IF(
metsol == 2)
THEN
10607 ELSE IF(
metsol == 3)
THEN
10610 ELSE IF(
metsol == 4)
THEN
10613 ELSE IF(
metsol == 5)
THEN
10616 ELSE IF(
metsol == 6)
THEN
10628 WRITE(
lunlog,*)
'Checking feasibility of parameters:'
10629 WRITE(*,*)
'Checking feasibility of parameters:'
10630 CALL feasib(concut,iact)
10632 WRITE(*,102) concut
10633 WRITE(*,*)
' parameters are made feasible'
10634 WRITE(
lunlog,102) concut
10635 WRITE(
lunlog,*)
' parameters are made feasible'
10637 WRITE(*,*)
' parameters are feasible (i.e. satisfy constraints)'
10638 WRITE(
lunlog,*)
' parameters are feasible (i.e. satisfy constraints)'
10646 WRITE(*,*)
'Reading files and accumulating vectors/matrices ...'
10650 WRITE(
lunlog,*)
'Reading files and accumulating vectors/matrices ...'
10667 IF(jcalcm+1 /= 0)
THEN
10676 IF(
iterat /= litera)
THEN
10697 IF (iabs(jcalcm) <= 1)
THEN
10710 IF(3*nrej >
nrecal)
THEN
10712 WRITE(*,*)
'Data records rejected in previous loop: '
10714 WRITE(*,*)
'Too many rejects (>33.3%) - stop'
10715 CALL peend(26,
'Aborted, too many rejects')
10722 IF(abs(
icalcm) == 1)
THEN
10734 ELSE IF(
metsol == 2)
THEN
10736 ELSE IF(
metsol == 3)
THEN
10738 ELSE IF(
metsol == 4)
THEN
10740 ELSE IF(
metsol == 5)
THEN
10742 ELSE IF(
metsol == 6)
THEN
10743 WRITE(*,*)
'... reserved for GMRES (not yet!)'
10746 ELSE IF(
metsol == 7)
THEN
10748 ELSE IF(
metsol == 8)
THEN
10751 ELSE IF(
metsol == 9)
THEN
10765 CALL feasib(concut,iact)
10777 angras=real(db/sqrt(db1*db2),mps)
10778 dbsig=16.0_mpd*sqrt(max(db1,db2))*epsilon(db)
10784 lsflag=lsflag .AND. (db > dbsig)
10791 ELSE IF(
metsol == 2)
THEN
10794 ELSE IF(
metsol == 3)
THEN
10797 ELSE IF(
metsol == 4)
THEN
10800 ELSE IF(
metsol == 5)
THEN
10803 ELSE IF(
metsol == 6)
THEN
10813 IF(db <= -dbsig)
THEN
10814 WRITE(*,*)
'Function not decreasing:',db
10815 IF(db > -1.0e-3_mpd)
THEN
10817 IF (iagain <= 1)
THEN
10818 WRITE(*,*)
'... again matrix calculation'
10822 WRITE(*,*)
'... aborting iterations'
10826 WRITE(*,*)
'... stopping iterations'
10849 IF (
nloopn == nloopsol)
THEN
10855 stepl=real(stp,mps)
10865 WRITE(*,*)
'Result vector containes ', nan,
' NaNs - stop'
10866 CALL peend(25,
'Aborted, result vector contains NaNs')
10873 WRITE(*,*)
'Subito! Exit after first step.'
10878 WRITE(*,*)
'INFO=0 should not happen (line search input err)'
10879 IF (iagain <= 0)
THEN
10884 IF(info < 0 .OR.
nloopn == nloopsol) cycle
10888 CALL feasib(concut,iact)
10889 IF(iact /= 0.OR.
chicut > 1.0)
THEN
10904 IF(sum(
nrejec) /= 0)
THEN
10906 WRITE(*,*)
'Data records rejected in last loop: '
10928 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
10929 IF (nfit > npar)
THEN
10932 WRITE(
lunlog,*)
'Inverse of global matrix from LDLt factorization'
10937#ifdef SCOREP_USER_ENABLE
10938 scorep_user_region_by_name_begin(
"UR_dsptri", scorep_user_region_type_common)
10941 IF(infolp /= 0) print *,
' DSPTRI failed: ', infolp
10942#ifdef SCOREP_USER_ENABLE
10943 scorep_user_region_by_name_end(
"UR_dsptri")
10949#ifdef SCOREP_USER_ENABLE
10950 scorep_user_region_by_name_begin(
"UR_dsytri", scorep_user_region_type_common)
10952 CALL dsytri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(nfit,mpl),&
10954 IF(infolp /= 0) print *,
' DSYTRI failed: ', infolp
10955#ifdef SCOREP_USER_ENABLE
10956 scorep_user_region_by_name_end(
"UR_dsytri")
10963 WRITE(
lunlog,*)
'Inverse of global matrix from LLt factorization'
10968#ifdef SCOREP_USER_ENABLE
10969 scorep_user_region_by_name_begin(
"UR_dpptri", scorep_user_region_type_common)
10971 CALL dpptri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),infolp)
10972 IF(infolp /= 0) print *,
' DPPTRI failed: ', infolp
10973#ifdef SCOREP_USER_ENABLE
10974 scorep_user_region_by_name_end(
"UR_dpptri")
10979#ifdef SCOREP_USER_ENABLE
10980 scorep_user_region_by_name_begin(
"UR_dpotri", scorep_user_region_type_common)
10982 CALL dpotri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(npar,mpl),infolp)
10983 IF(infolp /= 0) print *,
' DPOTRI failed: ', infolp
10984#ifdef SCOREP_USER_ENABLE
10985 scorep_user_region_by_name_end(
"UR_dpotri")
11003 DO i=npar-ncon+1,npar
11010 WRITE(
lunlog,*)
'Expansion of global matrix (A->Q*A*Q^t)'
11026 catio=real(dratio,mps)
11030 mrati=nint(100.0*catio,mpi)
11034 IF (
nfilw <= 0)
THEN
11035 WRITE(lunp,*)
'Sum(Chi^2)/Sum(Ndf) =',
fvalue
11037 WRITE(lunp,*)
' =',dratio
11039 WRITE(lunp,*)
'Sum(W*Chi^2)/Sum(Ndf)/<W> =',
fvalue
11041 WRITE(lunp,*)
' /',dwmean
11042 WRITE(lunp,*)
' =',dratio
11046 ' with correction for down-weighting ',catio
11067 nrati=nint(10000.0*real(nrej,mps)/real(
nrecal,mps),mpi)
11068 WRITE(crjrat,197) 0.01_mpd*real(nrati,mpd)
11069 nfaci=nint(100.0*sqrt(catio),mpi)
11071 WRITE(cratio,197) 0.01_mpd*real(mrati,mpd)
11072 WRITE(cfacin,197) 0.01_mpd*real(nfaci,mpd)
11075 IF(mrati < 90.OR.mrati > 110) warner=.true.
11076 IF(nrati > 100) warner=.true.
11077 IF(
ncgbe /= 0) warner=.true.
11079 IF(
nalow /= 0) warners=.true.
11081 IF(
nmiss1 /= 0) warnerss=.true.
11082 IF(iagain /= 0) warnerss=.true.
11083 IF(
ndefec /= 0) warnerss=.true.
11084 IF(
ndefpg /= 0) warnerss=.true.
11086 IF(
nrderr /= 0) warners3=.true.
11088 IF(warner.OR.warners.OR.warnerss.Or.warners3)
THEN
11091 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
11092 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
11093 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
11094 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
11095 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
11096 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
11097 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
11099 IF(mrati < 90.OR.mrati > 110)
THEN
11101 WRITE(*,*)
' Chi^2/Ndf = ',cratio,
' (should be close to 1)'
11102 WRITE(*,*)
' => multiply all input standard ', &
11103 'deviations by factor',cfacin
11106 IF(nrati > 100)
THEN
11108 WRITE(*,*)
' Fraction of rejects =',crjrat,
' %', &
11109 ' (should be far below 1 %)'
11110 WRITE(*,*)
' => please provide correct mille data'
11114 IF(iagain /= 0)
THEN
11116 WRITE(*,*)
' Matrix not positiv definite '// &
11117 '(function not decreasing)'
11118 WRITE(*,*)
' => please provide correct mille data'
11123 WRITE(*,*)
' Rank defect =',
ndefec, &
11124 ' for global matrix, should be 0'
11125 WRITE(*,*)
' => please provide correct mille data'
11130 WRITE(*,*)
' Rank defect for',
ndefpg, &
11131 ' parameter groups, should be 0'
11132 WRITE(*,*)
' => please provide correct mille data'
11137 WRITE(*,*)
' Rank defect =',
nmiss1, &
11138 ' for constraint equations, should be 0'
11139 WRITE(*,*)
' => please correct constraint definition'
11142 IF(
ncgbe /= 0)
THEN
11144 WRITE(*,*)
' Number of empty constraints =',
ncgbe,
', should be 0'
11145 WRITE(*,*)
' => please check constraint definition, mille data'
11148 IF(
nxlow /= 0)
THEN
11150 WRITE(*,*)
' Possible rank defects =',
nxlow,
' for global matrix'
11151 WRITE(*,*)
' (too few accepted entries)'
11152 WRITE(*,*)
' => please check mille data and ENTRIES cut'
11155 IF(
nalow /= 0)
THEN
11157 WRITE(*,*)
' Possible bad elements =',
nalow,
' in global vector'
11158 WRITE(*,*)
' (toos few accepted entries)'
11159 IF(
ipcntr > 0)
WRITE(*,*)
' (indicated in millepede.res by counts<0)'
11160 WRITE(*,*)
' => please check mille data and ENTRIES cut'
11165 WRITE(*,*)
' Binary file(s) with read errors =',
nrderr,
' (treated as EOF)'
11166 WRITE(*,*)
' => please check mille data'
11170 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
11171 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
11172 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
11173 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
11174 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
11175 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
11176 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
11187 ELSE IF(
metsol == 2)
THEN
11189 ELSE IF(
metsol == 3)
THEN
11195 IF(labelg == 0) cycle
11196 itgbi=inone(labelg)
11199 IF(ivgbi < 0) ivgbi=0
11200 IF(ivgbi == 0) cycle
11209 ELSE IF(
metsol == 6)
THEN
11212 ELSE IF(
metsol == 7)
THEN
11220 CALL peend(4,
'Ended with severe warnings (bad binary file(s))')
11221 ELSE IF (warnerss)
THEN
11222 CALL peend(3,
'Ended with severe warnings (bad global matrix)')
11223 ELSE IF (warners)
THEN
11224 CALL peend(2,
'Ended with severe warnings (insufficient measurements)')
11225 ELSE IF (warner)
THEN
11226 CALL peend(1,
'Ended with warnings (bad measurements)')
11228 CALL peend(0,
'Ended normally')
11231102
FORMAT(
' Call FEASIB with cut=',g10.3)
11249 INTEGER(mpi) :: kfl
11250 INTEGER(mpi) :: kmin
11251 INTEGER(mpi) :: kmax
11252 INTEGER(mpi) :: nrc
11253 INTEGER(mpl) :: nrej
11259 REAL(mpd) :: sumallw
11260 REAL(mpd) :: sumrejw
11262 sumallw=0.; sumrejw=0.;
11271 sumallw=sumallw+real(nrc,mpd)*
wfd(kfl)
11272 sumrejw=sumrejw+real(nrej,mpd)*
wfd(kfl)
11273 frac=real(nrej,mps)/real(nrc,mps)
11274 IF (frac > fmax)
THEN
11278 IF (frac < fmin)
THEN
11285 WRITE(*,
"(' Weighted fraction =',F8.2,' %')") 100.*sumrejw/sumallw
11286 IF (
nfilb > 1)
THEN
11287 WRITE(*,
"(' File with max. fraction ',I6,' :',F8.2,' %')") kmax, 100.*fmax
11288 WRITE(*,
"(' File with min. fraction ',I6,' :',F8.2,' %')") kmin, 100.*fmin
11314 INTEGER(mpi) :: iargc
11317 INTEGER(mpi) :: ierrf
11318 INTEGER(mpi) :: ieq
11319 INTEGER(mpi) :: ifilb
11320 INTEGER(mpi) :: ioff
11321 INTEGER(mpi) :: iopt
11322 INTEGER(mpi) :: ios
11323 INTEGER(mpi) :: iosum
11326 INTEGER(mpi) :: mat
11327 INTEGER(mpi) :: nab
11328 INTEGER(mpi) :: nline
11329 INTEGER(mpi) :: npat
11330 INTEGER(mpi) :: ntext
11332 INTEGER(mpi) :: nuf
11333 INTEGER(mpi) :: nums
11334 INTEGER(mpi) :: nufile
11335 INTEGER(mpi) :: lenfileInfo
11336 INTEGER(mpi) :: lenFileNames
11337 INTEGER(mpi) :: matint
11338 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: vecfileInfo
11339 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: tempArray
11340 INTEGER(mpl) :: rows
11341 INTEGER(mpl) :: cols
11342 INTEGER(mpl) :: newcols
11343 INTEGER(mpl) :: length
11345 CHARACTER (LEN=1024) :: text
11346 CHARACTER (LEN=1024) :: fname
11347 CHARACTER (LEN=14) :: bite(3)
11348 CHARACTER (LEN=32) :: keystx
11349 INTEGER(mpi),
PARAMETER :: mnum=100
11350 REAL(mpd) :: dnum(mnum)
11354 SUBROUTINE initc(nfiles)
BIND(c)
11356 INTEGER(c_int),
INTENT(IN),
VALUE :: nfiles
11357 END SUBROUTINE initc
11362 DATA bite/
'C_binary',
'text ',
'Fortran_binary'/
11377 WRITE(*,*)
'Command line options: '
11378 WRITE(*,*)
'--------------------- '
11380 CALL getarg(i,text)
11381 CALL rltext(text,ia,ib,nab)
11382 WRITE(*,101) i,text(1:nab)
11383 IF(text(ia:ia) /=
'-')
THEN
11384 nu=nufile(text(ia:ib))
11387 WRITE(*,*)
'Second text file in command line - stop'
11388 CALL peend(12,
'Aborted, second text file in command line')
11394 WRITE(*,*)
'Open error for file:',text(ia:ib),
' - stop'
11395 CALL peend(16,
'Aborted, open error for file')
11396 IF(text(ia:ia) /=
'/')
THEN
11397 CALL getenv(
'PWD',text)
11398 CALL rltext(text,ia,ib,nab)
11399 WRITE(*,*)
'PWD:',text(ia:ib)
11404 IF(index(text(ia:ib),
'b') /= 0)
THEN
11406 WRITE(*,*)
'Debugging requested'
11408 it=index(text(ia:ib),
't')
11411 ieq=index(text(ia+it:ib),
'=')+it
11412 IF (it /= ieq)
THEN
11413 IF (index(text(ia+ieq:ib),
'SL0' ) /= 0)
ictest=2
11414 IF (index(text(ia+ieq:ib),
'SLE' ) /= 0)
ictest=3
11415 IF (index(text(ia+ieq:ib),
'BP' ) /= 0)
ictest=4
11416 IF (index(text(ia+ieq:ib),
'BRLF') /= 0)
ictest=5
11417 IF (index(text(ia+ieq:ib),
'BRLC') /= 0)
ictest=6
11420 IF(index(text(ia:ib),
's') /= 0)
isubit=1
11421 IF(index(text(ia:ib),
'f') /= 0)
iforce=1
11422 IF(index(text(ia:ib),
'c') /= 0)
icheck=1
11423 IF(index(text(ia:ib),
'C') /= 0)
icheck=2
11425 IF(i == iargc())
WRITE(*,*)
'--------------------- '
11446 CALL rltext(text,ia,ib,nab)
11447 nu=nufile(text(ia:ib))
11451 CALL peend(10,
'Aborted, no steering file')
11452 stop
'in FILETC: no steering file. .'
11465 WRITE(*,*)
'Listing of steering file: ',
filnam(1:
nfnam)
11466 WRITE(*,*)
'-------------------------'
11469 WRITE(*,*)
'Open error for steering file - stop'
11470 CALL peend(11,
'Aborted, open error for steering file')
11471 IF(
filnam(1:1) /=
'/')
THEN
11472 CALL getenv(
'PWD',text)
11473 CALL rltext(text,ia,ib,nab)
11474 WRITE(*,*)
'PWD:',text(ia:ib)
11483 rows=6; cols=lenfileinfo
11484 CALL mpalloc(vecfileinfo,rows,cols,
'file info from steering')
11487 READ(10,102,iostat=ierrf) text
11488 IF (ierrf < 0)
EXIT
11489 CALL rltext(text,ia,ib,nab)
11491 IF(nline <= 50)
THEN
11492 WRITE(*,101) nline,text(1:nab)
11493 IF(nline == 50)
WRITE(*,*)
' ...'
11497 CALL rltext(text,ia,ib,nab)
11498 IF(ib == ia+2)
THEN
11499 mat=matint(text(ia:ib),
'end',npat,ntext)
11500 IF(mat == max(npat,ntext))
THEN
11503 WRITE(*,*)
' end-statement after',nline,
' text lines'
11508 keystx=
'fortranfiles'
11509 mat=matint(text(ia:ib),keystx,npat,ntext)
11510 IF(mat == max(npat,ntext))
THEN
11517 mat=matint(text(ia:ib),keystx,npat,ntext)
11518 IF(mat == max(npat,ntext))
THEN
11524 keystx=
'closeandreopen'
11525 mat=matint(text(ia:ib),keystx,npat,ntext)
11526 IF(mat == max(npat,ntext))
THEN
11534 iopt=index(text(ia:ib),
' -- ')
11535 IF (iopt > 0) ie=iopt-1
11538 nu=nufile(text(ia:ie))
11540 IF (
nfiles == lenfileinfo)
THEN
11541 CALL mpalloc(temparray,rows,cols,
'temp file info from steering')
11542 temparray=vecfileinfo
11544 lenfileinfo=lenfileinfo*2
11545 newcols=lenfileinfo
11546 CALL mpalloc(vecfileinfo,rows,newcols,
'file info from steering')
11547 vecfileinfo(:,1:cols)=temparray(:,1:cols)
11553 lenfilenames=lenfilenames+ie-ia+1
11554 vecfileinfo(1,
nfiles)=nline
11555 vecfileinfo(2,
nfiles)=nu
11556 vecfileinfo(3,
nfiles)=ia
11557 vecfileinfo(4,
nfiles)=ie
11558 vecfileinfo(5,
nfiles)=iopt
11559 vecfileinfo(6,
nfiles)=ib
11569 CALL mpalloc(
nfd,length,
'file line (in steering)')
11572 length=lenfilenames
11578 READ(10,102,iostat=ierrf) text
11579 IF (ierrf < 0)
EXIT
11581 IF (nline == vecfileinfo(1,i))
THEN
11582 nfd(i)=vecfileinfo(1,i)
11583 mfd(i)=vecfileinfo(2,i)
11584 ia=vecfileinfo(3,i)-1
11585 lfd(i)=vecfileinfo(4,i)-ia
11587 tfd(ioff+k)=text(ia+k:ia+k)
11592 IF (vecfileinfo(5,i) > 0)
THEN
11593 CALL ratext(text(vecfileinfo(5,i)+4:vecfileinfo(6,i)),nums,dnum,mnum)
11594 IF (nums > 0)
ofd(i)=real(dnum(1),mps)
11604 CALL mpalloc(
ifd,length,
'integrated record numbers (=offset)')
11605 CALL mpalloc(
jfd,length,
'number of accepted records')
11606 CALL mpalloc(
kfd,rows,length,
'number of records in file, file order')
11611 CALL mpalloc(
sfd,rows,length,
'start, end of file name in TFD')
11612 CALL mpalloc(
yfd,length,
'modification date')
11615 WRITE(*,*)
'-------------------------'
11621 WRITE(*,*)
'Table of files:'
11622 WRITE(*,*)
'---------------'
11625 WRITE(8,*)
'Text and data files:'
11629 fname(k:k)=
tfd(ioff+k)
11632 IF (
mprint > 1)
WRITE(*,103) i,bite(
mfd(i)),fname(1:
lfd(i))
11633 WRITE(8,103) i,bite(
mfd(i)),fname(1:
lfd(i))
11637 WRITE(*,*)
'---------------'
11651 IF(
mfd(i) == 3)
THEN
11675 IF(
mfd(i) == 1)
THEN
11696 WRITE(*,*)
'Opening of C-files not supported.'
11712 IF(iosum /= 0)
THEN
11713 CALL peend(15,
'Aborted, open error(s) for binary files')
11714 stop
'FILETC: open error '
11716 IF(
nfilb == 0)
THEN
11717 CALL peend(14,
'Aborted, no binary files')
11718 stop
'FILETC: no binary files '
11721 WRITE(*,*)
nfilb,
' binary files opened'
11723 WRITE(*,*)
nfilb,
' binary files opened and closed'
11727103
FORMAT(i3,2x,a14,3x,a)
11790 INTEGER(mpi) :: ierrf
11791 INTEGER(mpi) :: ioff
11792 INTEGER(mpi) :: ios
11793 INTEGER(mpi) :: iosum
11795 INTEGER(mpi) :: mat
11796 INTEGER(mpi) :: nab
11797 INTEGER(mpi) :: nfiln
11798 INTEGER(mpi) :: nline
11799 INTEGER(mpi) :: nlinmx
11800 INTEGER(mpi) :: npat
11801 INTEGER(mpi) :: ntext
11802 INTEGER(mpi) :: matint
11806 CHARACTER (LEN=1024) :: text
11807 CHARACTER (LEN=1024) :: fname
11810 WRITE(*,*)
'Processing text files ...'
11823 IF(
mfd(i) /= 2) cycle
11825 fname(k:k)=
tfd(ia+k)
11827 WRITE(*,*)
'File ',fname(1:
lfd(i))
11828 IF (
mprint > 1)
WRITE(*,*)
' '
11829 OPEN(10,file=fname(1:
lfd(i)),iostat=ios,form=
'FORMATTED')
11831 WRITE(*,*)
'Open error for file ',fname(1:
lfd(i))
11841 READ(10,102,iostat=ierrf) text
11842 IF (ierrf < 0)
THEN
11845 WRITE(*,*)
' end-of-file after',nline,
' text lines'
11849 IF(nline <= nlinmx.AND.
mprint > 1)
THEN
11850 CALL rltext(text,ia,ib,nab)
11852 WRITE(*,101) nline,text(1:nab)
11853 IF(nline == nlinmx)
WRITE(*,*)
' ...'
11856 CALL rltext(text,ia,ib,nab)
11857 IF(ib == ia+2)
THEN
11858 mat=matint(text(ia:ib),
'end',npat,ntext)
11859 IF(mat == max(npat,ntext))
THEN
11862 WRITE(*,*)
' end-statement after',nline,
' text lines'
11868 IF(nfiln <=
nfiles)
THEN
11869 IF(nline ==
nfd(nfiln))
THEN
11884 IF(iosum /= 0)
THEN
11885 CALL peend(16,
'Aborted, open error(s) for text files')
11886 stop
'FILETX: open error(s) in text files '
11889 WRITE(*,*)
'... end of text file processing.'
11894 WRITE(*,*)
lunkno,
' unknown keywords in steering files, ', &
11895 'or file non-existing,'
11896 WRITE(*,*)
' see above!'
11897 WRITE(*,*)
'------------> stop'
11899 CALL peend(13,
'Aborted, unknown keywords in steering file')
11908 ELSE IF(
matsto == 1)
THEN
11910 ELSE IF(
matsto == 2)
THEN
11913 ELSE IF(
metsol == 1)
THEN
11915 ELSE IF(
metsol == 2)
THEN
11917 ELSE IF(
metsol == 3)
THEN
11919 ELSE IF(
metsol == 4)
THEN
11921 ELSE IF(
metsol == 5)
THEN
11923 ELSE IF(
metsol == 6)
THEN
11926 ELSE IF(
metsol == 7)
THEN
11928 ELSE IF(
metsol == 8)
THEN
11931 ELSE IF(
metsol == 9)
THEN
11936 WRITE(*,*)
'MINRES forced with sparse matrix!'
11938 WRITE(*,*)
'MINRES forced with sparse matrix!'
11940 WRITE(*,*)
'MINRES forced with sparse matrix!'
11945 WRITE(*,*)
'MINRES forced with sparse matrix!'
11947 WRITE(*,*)
'MINRES forced with sparse matrix!'
11949 WRITE(*,*)
'MINRES forced with sparse matrix!'
11957 WRITE(*,*)
'Solution method and matrix-storage mode:'
11959 WRITE(*,*)
' METSOL = 1: matrix inversion'
11960 ELSE IF(
metsol == 2)
THEN
11961 WRITE(*,*)
' METSOL = 2: diagonalization'
11962 ELSE IF(
metsol == 3)
THEN
11963 WRITE(*,*)
' METSOL = 3: decomposition'
11964 ELSE IF(
metsol == 4)
THEN
11965 WRITE(*,*)
' METSOL = 4: MINRES'
11966 ELSE IF(
metsol == 5)
THEN
11967 WRITE(*,*)
' METSOL = 5: MINRES-QLP'
11968 ELSE IF(
metsol == 6)
THEN
11969 WRITE(*,*)
' METSOL = 6: GMRES (-> MINRES)'
11971 ELSE IF(
metsol == 7)
THEN
11972 WRITE(*,*)
' METSOL = 7: LAPACK factorization'
11973 ELSE IF(
metsol == 8)
THEN
11974 WRITE(*,*)
' METSOL = 8: LAPACK factorization'
11976 ELSE IF(
metsol == 9)
THEN
11977 WRITE(*,*)
' METSOL = 9: Intel oneMKL PARDISO'
11982 WRITE(*,*)
' with',
mitera,
' iterations'
11985 WRITE(*,*)
' MATSTO = 0: unpacked symmetric matrix, ',
'n*n elements'
11986 ELSEIF(
matsto == 1)
THEN
11987 WRITE(*,*)
' MATSTO = 1: full symmetric matrix, ',
'(n*n+n)/2 elements'
11988 ELSE IF(
matsto == 2)
THEN
11989 WRITE(*,*)
' MATSTO = 2: sparse matrix (custom)'
11990 ELSE IF(
matsto == 3)
THEN
11992 WRITE(*,*)
' MATSTO = 3: sparse matrix (upper triangle, CSR3)'
11994 WRITE(*,*)
' MATSTO = 3: sparse matrix (upper triangle, BSR3)'
11998 WRITE(*,*)
' and band matrix, width',
mbandw
12002 WRITE(*,*)
'Chi square cut equiv 3 st.dev applied ...'
12003 WRITE(*,*)
' in first iteration with factor',
chicut
12004 WRITE(*,*)
' in second iteration with factor',
chirem
12005 WRITE(*,*)
' (reduced by sqrt in next iterations)'
12009 WRITE(*,*)
' Down-weighting of outliers in',
lhuber,
' iterations'
12010 WRITE(*,*)
' Cut on downweight fraction',
dwcut
12013 WRITE(*,*)
'Iterations (solutions) with line search:'
12022 WRITE(*,*)
' All with Chi square cut scaling factor <= 1.'
12053 INTEGER(mpi) :: ios
12057 INTEGER(mpi) :: npat
12058 INTEGER(mpi) :: ntext
12059 INTEGER(mpi) :: nuprae
12062 CHARACTER (LEN=*),
INTENT(INOUT) :: fname
12068 IF(len(fname) > 5)
THEN
12069 IF(fname(1:5) ==
'rfio:') nuprae=1
12070 IF(fname(1:5) ==
'dcap:') nuprae=2
12071 IF(fname(1:5) ==
'root:') nuprae=3
12073 IF(nuprae == 0)
THEN
12074 INQUIRE(file=fname,iostat=ios,exist=ex)
12075 IF(ios /= 0)
nufile=-abs(ios)
12076 IF(ios /= 0)
RETURN
12077 ELSE IF(nuprae == 1)
THEN
12090 nm=
matint(
'xt',fname(l1:ll),npat,ntext)
12093 nm=
matint(
'tx',fname(l1:ll),npat,ntext)
12114 INTEGER(mpi) :: ier
12115 INTEGER(mpi) :: iomp
12118 INTEGER(mpi) :: kkey
12119 INTEGER(mpi) :: label
12120 INTEGER(mpi) :: lkey
12121 INTEGER(mpi) :: mat
12122 INTEGER(mpi) :: miter
12123 INTEGER(mpi) :: nab
12124 INTEGER(mpi) :: nkey
12125 INTEGER(mpi) :: nkeys
12127 INTEGER(mpi) :: nmeth
12128 INTEGER(mpi) :: npat
12129 INTEGER(mpi) :: ntext
12130 INTEGER(mpi) :: nums
12131 INTEGER(mpi) :: matint
12133 CHARACTER (LEN=*),
INTENT(IN) :: text
12134 INTEGER(mpi),
INTENT(IN) :: nline
12138 parameter(nkeys=7,nmeth=10)
12140 parameter(nkeys=6,nmeth=9)
12143 parameter(nkeys=6,nmeth=7)
12145 CHARACTER (LEN=16) :: methxt(nmeth)
12146 CHARACTER (LEN=16) :: keylst(nkeys)
12147 CHARACTER (LEN=32) :: keywrd
12148 CHARACTER (LEN=32) :: keystx
12149 CHARACTER (LEN=itemCLen) :: ctext
12150 INTEGER(mpi),
PARAMETER :: mnum=100
12151 REAL(mpd) :: dnum(mnum)
12154 INTEGER(mpi) :: ipvs
12157 INTEGER(mpi) :: lpvs
12161 SUBROUTINE additem(length,list,label,value)
12163 INTEGER(mpi),
INTENT(IN OUT) :: length
12164 TYPE(listitem),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12165 INTEGER(mpi),
INTENT(IN) :: label
12166 REAL(mpd),
INTENT(IN) :: value
12168 SUBROUTINE additemc(length,list,label,text)
12170 INTEGER(mpi),
INTENT(IN OUT) :: length
12171 TYPE(listitemc),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12172 INTEGER(mpi),
INTENT(IN) :: label
12173 CHARACTER(LEN = itemCLen),
INTENT(IN) :: text
12175 SUBROUTINE additemi(length,list,label,ivalue)
12177 INTEGER(mpi),
INTENT(IN OUT) :: length
12178 TYPE(listitemi),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12179 INTEGER(mpi),
INTENT(IN) :: label
12180 INTEGER(mpi),
INTENT(IN) :: ivalue
12187 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment',
'pardiso'/
12188 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
12189 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition',
'fullLAPACK',
'unpackedLAPACK', &
12192 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment'/
12193 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
12194 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition',
'fullLAPACK',
'unpackedLAPACK'/
12197 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment'/
12198 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
12199 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition'/
12205 CALL rltext(text,ia,ib,nab)
12206 IF(nab == 0)
GOTO 10
12207 CALL ratext(text(1:nab),nums,dnum,mnum)
12209 IF(nums /= 0) nkey=0
12217 keystx=keylst(nkey)
12218 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12219 IF(100*mat >= 80*max(npat,ntext))
GO TO 10
12225 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12226 IF(100*mat >= 80*max(npat,ntext))
THEN
12228 IF(nums > 0) mprint=nint(dnum(1),mpi)
12233 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12234 IF(100*mat >= 80*max(npat,ntext))
THEN
12237 IF(nums > 0) mdebug=nint(dnum(1),mpi)
12238 IF(nums > 1) mdebg2=nint(dnum(2),mpi)
12243 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12244 IF(100*mat >= 80*max(npat,ntext))
THEN
12245 IF(nums > 0 .AND. dnum(1) > 0.5) mreqenf=nint(dnum(1),mpi)
12246 IF(nums > 1 .AND. dnum(2) > 0.5) mreqena=nint(dnum(2),mpi)
12247 IF(nums > 2 .AND. dnum(3) > 0.5) iteren=nint(dnum(1)*dnum(3),mpi)
12251 keystx=
'printrecord'
12252 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12253 IF(100*mat >= 80*max(npat,ntext))
THEN
12254 IF(nums > 0) nrecpr=nint(dnum(1),mpi)
12255 IF(nums > 1) nrecp2=nint(dnum(2),mpi)
12260 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12261 IF(100*mat >= 80*max(npat,ntext))
THEN
12262 IF (nums > 0.AND.dnum(1) > 0.) mxrec=nint(dnum(1),mpi)
12267 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12268 IF(100*mat >= 80*max(npat,ntext))
THEN
12269 IF (nums > 0.AND.dnum(1) >= 0.) ncache=nint(dnum(1),mpi)
12270 IF (nums == 2.AND.dnum(2) > 0..AND.dnum(2) <= 1.0) &
12271 fcache(1)=real(dnum(2),mps)
12272 IF (nums >= 4)
THEN
12274 fcache(k)=real(dnum(k+1),mps)
12281 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12282 IF(100*mat >= 80*max(npat,ntext))
THEN
12287 chicut=real(dnum(1),mps)
12288 IF(chicut < 1.0) chicut=-1.0
12292 chirem=real(dnum(2),mps)
12293 IF(chirem < 1.0) chirem=1.0
12294 IF(chicut >= 1.0) chirem=min(chirem,chicut)
12302 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12303 IF(100*mat >= 80*max(npat,ntext))
THEN
12304 IF(nums > 0) chhuge=real(dnum(1),mps)
12305 IF(chhuge < 1.0) chhuge=1.0
12310 keystx=
'linesearch'
12311 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12312 IF(100*mat >= 80*max(npat,ntext))
THEN
12313 IF(nums > 0) lsearch=nint(dnum(1),mpi)
12318 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12319 IF(100*mat >= 80*max(npat,ntext))
THEN
12320 IF(nums > 0) lfitnp=nint(dnum(1),mpi)
12321 IF(nums > 1) lfitbb=nint(dnum(2),mpi)
12325 keystx=
'regularization'
12326 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12327 IF(100*mat >= 80*max(npat,ntext))
THEN
12329 regula=real(dnum(1),mps)
12330 IF(nums >= 2) regpre=real(dnum(2),mps)
12334 keystx=
'regularisation'
12335 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12336 IF(100*mat >= 80*max(npat,ntext))
THEN
12338 regula=real(dnum(1),mps)
12339 IF(nums >= 2) regpre=real(dnum(2),mps)
12344 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12345 IF(100*mat >= 80*max(npat,ntext))
THEN
12346 regpre=real(dnum(1),mps)
12351 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12352 IF(100*mat >= 80*max(npat,ntext))
THEN
12353 matrit=nint(dnum(1),mpi)
12358 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12359 IF(100*mat >= 80*max(npat,ntext))
THEN
12361 IF (nums > 0.AND.dnum(1) > 0.) matmon=nint(dnum(1),mpi)
12366 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12367 IF(100*mat >= 80*max(npat,ntext))
THEN
12368 IF(nums > 0) mbandw=nint(dnum(1),mpi)
12369 IF(mbandw < 0) mbandw=-1
12370 IF(nums > 1) lprecm=nint(dnum(2),mpi)
12394 keystx=
'outlierdownweighting'
12395 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12396 IF(100*mat >= 80*max(npat,ntext))
THEN
12397 lhuber=nint(dnum(1),mpi)
12398 IF(lhuber > 0.AND.lhuber <= 2) lhuber=2
12402 keystx=
'dwfractioncut'
12403 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12404 IF(100*mat >= 80*max(npat,ntext))
THEN
12405 dwcut=real(dnum(1),mps)
12406 IF(dwcut > 0.5) dwcut=0.5
12410 keystx=
'maxlocalcond'
12411 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12412 IF(100*mat >= 80*max(npat,ntext))
THEN
12413 IF (nums > 0.AND.dnum(1) > 0.0) cndlmx=real(dnum(1),mps)
12418 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12419 IF(100*mat >= 80*max(npat,ntext))
THEN
12420 prange=abs(real(dnum(1),mps))
12425 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12426 IF(100*mat >= 80*max(npat,ntext))
THEN
12432 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12433 IF(100*mat >= 80*max(npat,ntext))
THEN
12438 keystx=
'memorydebug'
12439 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12440 IF(100*mat >= 80*max(npat,ntext))
THEN
12442 IF (nums > 0.AND.dnum(1) > 0.0) memdbg=nint(dnum(1),mpi)
12446 keystx=
'globalcorr'
12447 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12448 IF(100*mat >= 80*max(npat,ntext))
THEN
12453 keystx=
'printcounts'
12454 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12455 IF(100*mat >= 80*max(npat,ntext))
THEN
12457 IF (nums > 0) ipcntr=nint(dnum(1),mpi)
12461 keystx=
'weightedcons'
12462 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12463 IF(100*mat >= 80*max(npat,ntext))
THEN
12465 IF (nums > 0) iwcons=nint(dnum(1),mpi)
12469 keystx=
'skipemptycons'
12470 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12471 IF(100*mat >= 80*max(npat,ntext))
THEN
12476 keystx=
'resolveredundancycons'
12477 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12478 IF(100*mat >= 80*max(npat,ntext))
THEN
12483 keystx=
'withelimination'
12484 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12485 IF(100*mat >= 80*max(npat,ntext))
THEN
12490 keystx=
'postprocessing'
12491 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12492 IF(100*mat >= 80*max(npat,ntext))
THEN
12493 lenpostproc=ib-
keyb-1
12494 cpostproc(1:lenpostproc)=text(
keyb+2:ib)
12499 keystx=
'withLAPACKelimination'
12500 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12501 IF(100*mat >= 80*max(npat,ntext))
THEN
12507 keystx=
'withmultipliers'
12508 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12509 IF(100*mat >= 80*max(npat,ntext))
THEN
12514 keystx=
'checkinput'
12515 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12516 IF(100*mat >= 80*max(npat,ntext))
THEN
12518 IF (nums > 0) icheck=nint(dnum(1),mpi)
12522 keystx=
'checkparametergroups'
12523 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12524 IF(100*mat >= 80*max(npat,ntext))
THEN
12529 keystx=
'monitorresiduals'
12530 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12531 IF(100*mat >= 80*max(npat,ntext))
THEN
12533 IF (nums > 0) imonit=nint(dnum(1),mpi)
12534 IF (nums > 1) measbins=max(measbins,nint(dnum(2),mpi))
12538 keystx=
'monitorpulls'
12539 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12540 IF(100*mat >= 80*max(npat,ntext))
THEN
12543 IF (nums > 0) imonit=nint(dnum(1),mpi)
12544 IF (nums > 1) measbins=max(measbins,nint(dnum(2),mpi))
12548 keystx=
'monitorprogress'
12549 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12550 IF(100*mat >= 80*max(npat,ntext))
THEN
12553 IF (nums > 0) monpg1=max(1,nint(dnum(1),mpi))
12554 IF (nums > 1) monpg2=max(1,nint(dnum(2),mpi))
12558 keystx=
'scaleerrors'
12559 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12560 IF(100*mat >= 80*max(npat,ntext))
THEN
12562 IF (nums > 0) dscerr(1:2)=dnum(1)
12563 IF (nums > 1) dscerr(2)=dnum(2)
12567 keystx=
'iterateentries'
12568 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12569 IF(100*mat >= 80*max(npat,ntext))
THEN
12570 iteren=huge(iteren)
12571 IF (nums > 0) iteren=nint(dnum(1),mpi)
12576 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12577 IF(100*mat >= 80*max(npat,ntext))
THEN
12585 WRITE(*,*)
'WARNING: multithreading not available'
12591 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12592 IF(100*mat >= 80*max(npat,ntext))
THEN
12593 WRITE(*,*)
'WARNING: keyword COMPRESS is obsolete (compression is default)'
12605 keystx=
'countrecords'
12606 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12607 IF(100*mat >= 80*max(npat,ntext))
THEN
12613 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12614 IF(100*mat >= 80*max(npat,ntext).AND.mnrsel < 100)
THEN
12615 nl=min(nums,100-mnrsel)
12617 lbmnrs(mnrsel+k)=nint(dnum(k),mpi)
12623 keystx=
'pairentries'
12624 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12625 IF(100*mat >= 80*max(npat,ntext))
THEN
12628 IF (nums > 0.AND.dnum(1) > 0.0)
THEN
12629 mreqpe=nint(dnum(1),mpi)
12630 IF (nums >= 2.AND.dnum(2) >= dnum(1)) mhispe=nint(dnum(2),mpi)
12631 IF (nums >= 3.AND.dnum(3) >= dnum(1)) msngpe=nint(dnum(3),mpi)
12637 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12638 IF(100*mat >= 80*max(npat,ntext))
THEN
12639 wolfc1=real(dnum(1),mps)
12640 wolfc2=real(dnum(2),mps)
12647 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12648 IF(100*mat >= 80*max(npat,ntext))
THEN
12650 IF (dnum(1) < 1.0e-10_mpd.OR.dnum(1) > 1.0e-04_mpd)
THEN
12651 WRITE(*,*)
'ERROR: need 1.0D-10 <= MRESTL ', &
12652 '<= 1.0D-04, but get ', dnum(1)
12661 keystx=
'mrestranscond'
12662 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12663 IF(100*mat >= 80*max(npat,ntext))
THEN
12671 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12672 IF(100*mat >= 80*max(npat,ntext))
THEN
12674 mrmode = int(dnum(1),mpi)
12679 keystx=
'nofeasiblestart'
12680 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12681 IF(100*mat >= 80*max(npat,ntext))
THEN
12687 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12688 IF(100*mat >= 80*max(npat,ntext))
THEN
12693 keystx=
'readerroraseof'
12694 mat=matint(text(ia:ib),keystx,npat,ntext)
12695 IF(100*mat >= 80*max(npat,ntext))
THEN
12701 keystx=
'LAPACKwitherrors'
12702 mat=matint(text(ia:ib),keystx,npat,ntext)
12703 IF(100*mat >= 80*max(npat,ntext))
THEN
12708 keystx=
'debugPARDISO'
12709 mat=matint(text(ia:ib),keystx,npat,ntext)
12710 IF(100*mat >= 80*max(npat,ntext))
THEN
12715 keystx=
'blocksizePARDISO'
12716 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12717 IF(100*mat >= 80*max(npat,ntext).AND.mnrsel < 100)
THEN
12718 nl=min(nums,10-mpdbsz)
12720 IF (nint(dnum(k),mpi) > 0)
THEN
12721 IF (mpdbsz == 0)
THEN
12723 ipdbsz(mpdbsz)=nint(dnum(k),mpi)
12724 ELSE IF (nint(dnum(k),mpi) > ipdbsz(mpdbsz))
THEN
12726 ipdbsz(mpdbsz)=nint(dnum(k),mpi)
12734 keystx=
'fortranfiles'
12735 mat=matint(text(ia:ib),keystx,npat,ntext)
12736 IF(mat == max(npat,ntext))
RETURN
12739 mat=matint(text(ia:ib),keystx,npat,ntext)
12740 IF(mat == max(npat,ntext))
RETURN
12742 keystx=
'closeandreopen'
12743 mat=matint(text(ia:ib),keystx,npat,ntext)
12744 IF(mat == max(npat,ntext))
RETURN
12748 IF(nums /= 0) nkey=0
12751 WRITE(*,*)
'**************************************************'
12753 WRITE(*,*)
'Unknown keyword(s): ',text(1:min(nab,50))
12755 WRITE(*,*)
'**************************************************'
1277710
IF(nkey > 0)
THEN
12781 lpvs=nint(dnum(1),mpi)
12783 CALL additem(lenparameters,listparameters,lpvs,dnum(2))
12784 CALL additem(lenpresigmas,listpresigmas,lpvs,dnum(3))
12786 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12788 ELSE IF(nums /= 0)
THEN
12790 WRITE(*,*)
'Wrong text in line',nline
12791 WRITE(*,*)
'Status: new parameter'
12792 WRITE(*,*)
'> ',text(1:nab)
12794 ELSE IF(lkey == 3)
THEN
12796 IF(nums >= 1.AND.nums <= 2)
THEN
12798 CALL additem(lenconstraints,listconstraints,lpvs,dnum(1))
12800 IF(iwcons > 0) lpvs=-2
12802 IF(nums == 2) plvs=dnum(2)
12803 CALL additem(lenconstraints,listconstraints,lpvs,plvs)
12806 WRITE(*,*)
'Wrong text in line',nline
12807 WRITE(*,*)
'Status: new keyword constraint'
12808 WRITE(*,*)
'> ',text(1:nab)
12810 ELSE IF(lkey == 4)
THEN
12812 nummeasurements=nummeasurements+1
12814 CALL additem(lenmeasurements,listmeasurements,lpvs,dnum(1))
12816 CALL additem(lenmeasurements,listmeasurements,lpvs,dnum(2))
12819 WRITE(*,*)
'Wrong text in line',nline
12820 WRITE(*,*)
'Status: new keyword measurement'
12821 WRITE(*,*)
'> ',text(1:nab)
12823 ELSE IF(lkey == 5.AND.
keyb <
keyc)
THEN
12825 IF(nums >= 1) miter=nint(dnum(1),mpi)
12826 IF(miter >= 1) mitera=miter
12827 dflim=real(dnum(2),mps)
12831 mat=matint(text(
keyb+1:
keyc),keystx,npat,ntext)
12832 IF(100*mat >= 80*max(npat,ntext))
THEN
12836 ELSE IF(i == 2)
THEN
12839 ELSE IF(i == 3)
THEN
12842 ELSE IF(i == 4)
THEN
12845 ELSE IF(i == 5)
THEN
12848 ELSE IF(i == 6)
THEN
12851 ELSE IF(i == 7)
THEN
12855 ELSE IF(i == 8)
THEN
12858 ELSE IF(i == 9)
THEN
12862 ELSE IF(i == 10)
THEN
12871 ELSE IF(nkey == 0)
THEN
12874 lpvs=nint(dnum(1),mpi)
12876 CALL additem(lenparameters,listparameters,lpvs,dnum(2))
12877 CALL additem(lenpresigmas,listpresigmas,lpvs,dnum(3))
12879 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12881 ELSE IF(nums > 1.AND.nums < 3)
THEN
12883 WRITE(*,*)
'Wrong text in line',nline
12884 WRITE(*,*)
'Status continuation parameter'
12885 WRITE(*,*)
'> ',text(1:nab)
12888 ELSE IF(lkey == 3)
THEN
12891 label=nint(dnum(i),mpi)
12892 IF(label <= 0) ier=1
12894 IF(mod(nums,2) /= 0) ier=1
12897 lpvs=nint(dnum(i),mpi)
12899 CALL additem(lenconstraints,listconstraints,lpvs,plvs)
12903 WRITE(*,*)
'Wrong text in line',nline
12904 WRITE(*,*)
'Status continuation constraint'
12905 WRITE(*,*)
'> ',text(1:nab)
12908 ELSE IF(lkey == 4)
THEN
12912 label=nint(dnum(i),mpi)
12913 IF(label <= 0) ier=1
12915 IF(mod(nums,2) /= 0) ier=1
12919 lpvs=nint(dnum(i),mpi)
12921 CALL additem(lenmeasurements,listmeasurements,lpvs,plvs)
12925 WRITE(*,*)
'Wrong text in line',nline
12926 WRITE(*,*)
'Status continuation measurement'
12927 WRITE(*,*)
'> ',text(1:nab)
12929 ELSE IF(lkey == 6)
THEN
12931 lpvs=nint(dnum(1),mpi)
12935 IF (text(j:j) ==
' ')
EXIT
12938 CALL additemc(lencomments,listcomments,lpvs,ctext)
12940 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12942 ELSE IF(nums /= 0)
THEN
12944 WRITE(*,*)
'Wrong text in line',nline
12945 WRITE(*,*)
'Status: continuation comment'
12946 WRITE(*,*)
'> ',text(1:nab)
12950 ELSE IF(lkey == 7)
THEN
12953 label=nint(dnum(i),mpi)
12954 IF(label <= 0.OR.label > 64) ier=1
12956 IF(mod(nums,2) /= 0) ier=1
12960 lpvs=nint(dnum(i),mpi)
12961 ipvs=nint(dnum(i+1),mpi)
12962 CALL additemi(lenpardiso,listpardiso,lpvs,ipvs)
12966 WRITE(*,*)
'Wrong text in line',nline
12967 WRITE(*,*)
'Status continuation measurement'
12968 WRITE(*,*)
'> ',text(1:nab)
12987 INTEGER(mpi),
INTENT(IN OUT) :: length
12988 TYPE(
listitem),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12989 INTEGER(mpi),
INTENT(IN) :: label
12990 REAL(mpd),
INTENT(IN) :: value
12992 INTEGER(mpl) :: newSize
12993 INTEGER(mpl) :: oldSize
12994 TYPE(
listitem),
DIMENSION(:),
ALLOCATABLE :: tempList
12996 IF (label > 0.AND.
value == 0.)
RETURN
12997 IF (length == 0 )
THEN
12999 CALL mpalloc(list,newsize,
' list ')
13001 oldsize=
size(list,kind=
mpl)
13002 IF (length >= oldsize)
THEN
13003 newsize = oldsize + oldsize/5 + 100
13004 CALL mpalloc(templist,oldsize,
' temp. list ')
13007 CALL mpalloc(list,newsize,
' list ')
13008 list(1:oldsize)=templist(1:oldsize)
13013 list(length)%label=label
13014 list(length)%value=
value
13029 INTEGER(mpi),
INTENT(IN OUT) :: length
13030 TYPE(
listitemc),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
13031 INTEGER(mpi),
INTENT(IN) :: label
13032 CHARACTER(len = itemCLen),
INTENT(IN) :: text
13034 INTEGER(mpl) :: newSize
13035 INTEGER(mpl) :: oldSize
13036 TYPE(
listitemc),
DIMENSION(:),
ALLOCATABLE :: tempList
13038 IF (label > 0.AND.text ==
'')
RETURN
13039 IF (length == 0 )
THEN
13041 CALL mpalloc(list,newsize,
' list ')
13043 oldsize=
size(list,kind=
mpl)
13044 IF (length >= oldsize)
THEN
13045 newsize = oldsize + oldsize/5 + 100
13046 CALL mpalloc(templist,oldsize,
' temp. list ')
13049 CALL mpalloc(list,newsize,
' list ')
13050 list(1:oldsize)=templist(1:oldsize)
13055 list(length)%label=label
13056 list(length)%text=text
13071 INTEGER(mpi),
INTENT(IN OUT) :: length
13072 TYPE(
listitemi),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
13073 INTEGER(mpi),
INTENT(IN) :: label
13074 INTEGER(mpi),
INTENT(IN) :: ivalue
13076 INTEGER(mpl) :: newSize
13077 INTEGER(mpl) :: oldSize
13078 TYPE(
listitemi),
DIMENSION(:),
ALLOCATABLE :: tempList
13080 IF (length == 0 )
THEN
13082 CALL mpalloc(list,newsize,
' list ')
13084 oldsize=
size(list,kind=
mpl)
13085 IF (length >= oldsize)
THEN
13086 newsize = oldsize + oldsize/5 + 100
13087 CALL mpalloc(templist,oldsize,
' temp. list ')
13090 CALL mpalloc(list,newsize,
' list ')
13091 list(1:oldsize)=templist(1:oldsize)
13096 list(length)%label=label
13097 list(length)%ivalue=ivalue
13111 CHARACTER (LEN=*),
INTENT(IN) :: text
13112 CHARACTER (LEN=16) :: textc
13121 textl(ka:kb)=text(1:l)
13125 textc=text(1:l)//
'-end'
13133 textl(ka:kb)=textc(1:l)
13160 INTEGER(mpi),
INTENT(IN) :: lun
13161 CHARACTER (LEN=*),
INTENT(IN) :: fname
13162 CHARACTER (LEN=33) :: nafile
13163 CHARACTER (LEN=33) :: nbfile
13169 CALL peend(17,
'Aborted, file name too long')
13170 stop
'File name too long '
13173 nafile(l+1:l+1)=
'~'
13175 INQUIRE(file=nafile(1:l),exist=ex)
13177 INQUIRE(file=nafile(1:l+1),exist=ex)
13179 CALL system(
'rm '//nafile)
13182 nafile(l+1:l+1)=
' '
13183 CALL system(
'mv '//nafile//nbfile)
13185 OPEN(unit=lun,file=fname)
13196 REAL,
DIMENSION(2) :: ta
13216 IF(ncount > 1)
THEN
13218 nsecd1=int(delta,
mpi)
13220 minut1=nsecd1/60-60*nhour1
13221 secnd1=delta-60*(minut1+60*nhour1)
13223 nsecd2=int(delta,
mpi)
13225 minut2=nsecd2/60-60*nhour2
13226 secnd2=delta-60*(minut2+60*nhour2)
13227 WRITE(*,101) nhour1,minut1,secnd1, nhour2,minut2,secnd2
13232101
FORMAT(i4,
' h',i3,
' min',f5.1,
' sec total',18x,
'elapsed', &
13233 i4,
' h',i3,
' min',f5.1,
' sec')
13247 INTEGER(mpi),
INTENT(IN) :: icode
13248 CHARACTER (LEN=*),
INTENT(IN) :: cmessage
13250 CALL mvopen(9,
'millepede.end')
13251 WRITE(9,101) icode, cmessage
13252101
FORMAT(1x,i4,3x,a)
13256END SUBROUTINE peend
13268 INTEGER(mpi),
INTENT(IN) :: kfile
13269 INTEGER(mpi),
INTENT(IN) :: ithr
13270 INTEGER(mpi),
INTENT(OUT) :: ierr
13272 INTEGER(mpi),
DIMENSION(13) :: ibuff
13273 INTEGER(mpi) :: ioff
13274 INTEGER(mpi) :: ios
13276 INTEGER(mpi) :: lfn
13277 INTEGER(mpi) :: lun
13278 INTEGER(mpi) :: moddate
13279 CHARACTER (LEN=1024) :: fname
13280 CHARACTER (LEN=7) :: cfile
13285 SUBROUTINE openc(filename, lfn, lun, ios)
BIND(c)
13287 CHARACTER(kind=c_char),
DIMENSION(*),
INTENT(IN) :: filename
13288 INTEGER(c_int),
INTENT(IN),
VALUE :: lfn
13289 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13290 INTEGER(c_int),
INTENT(INOUT) :: ios
13291 END SUBROUTINE openc
13303 fname(k:k)=
tfd(ioff+k)
13308 IF(kfile <=
nfilf)
THEN
13311 OPEN(lun,file=fname(1:lfn),iostat=ios, form=
'UNFORMATTED')
13312 print *,
' lun ', lun, ios
13316 CALL openc(fname(1:lfn),lfn,lun,ios)
13318 WRITE(*,*)
'Opening of C-files not supported.'
13325 WRITE(*,*)
'Open error for file ',fname(1:lfn), ios
13326 IF (moddate /= 0)
THEN
13327 WRITE(cfile,
'(I7)') kfile
13328 CALL peend(15,
'Aborted, open error(s) for binary file ' // cfile)
13329 stop
'PEREAD: open error'
13334 ios=stat(fname(1:lfn),ibuff)
13338 WRITE(*,*)
'STAT error for file ',fname(1:lfn), ios
13342 IF (moddate /= 0)
THEN
13343 IF (ibuff(10) /= moddate)
THEN
13344 WRITE(cfile,
'(I7)') kfile
13345 CALL peend(19,
'Aborted, binary file modified (date) ' // cfile)
13346 stop
'PEREAD: file modified'
13349 yfd(kfile)=ibuff(10)
13364 INTEGER(mpi),
INTENT(IN) :: kfile
13365 INTEGER(mpi),
INTENT(IN) :: ithr
13367 INTEGER(mpi) :: lun
13371 SUBROUTINE closec(lun)
BIND(c)
13373 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13380 IF(kfile <=
nfilf)
THEN
13399 INTEGER(mpi),
INTENT(IN) :: kfile
13401 INTEGER(mpi) :: lun
13405 SUBROUTINE resetc(lun)
BIND(c)
13407 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13413 IF (kfile <=
nfilf)
THEN
13432 INTEGER(mpi) :: ipgrp
13433 INTEGER(mpi) :: irank
13434 INTEGER(mpi) :: isize
13435 INTEGER(mpi) :: ivoff
13436 INTEGER(mpi) :: itgbi
13438 INTEGER(mpi) :: msize
13439 INTEGER(mpi),
PARAMETER :: mxsize = 1000
13441 INTEGER(mpl):: length
13443 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: auxVectorD
13444 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: auxVectorI
13445 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: resParGroup
13446 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: blockParGroup
13454 IF (isize <= mxsize)
THEN
13455 msize=max(msize,isize)
13457 print *,
' CKPGRP: par. group', ipgrp,
' not checked -- too large: ', isize
13460 IF (msize == 0)
RETURN
13463 length=int(msize,mpl)*(int(msize,mpl)+1)/2
13464 CALL mpalloc(blockpargroup,length,
'(matrix) block for parameter groups (D)')
13466 CALL mpalloc(respargroup,length,
'residuals for parameter groups (D)')
13467 CALL mpalloc(auxvectori,length,
'auxiliary array (I)')
13468 CALL mpalloc(auxvectord,length,
'auxiliary array (D)')
13472 print *,
' CKPGRP par. group first label size rank'
13475 IF (isize > mxsize) cycle
13482 blockpargroup(ij)=matij(ivoff+i,ivoff+j)
13486 CALL sqminv(blockpargroup,respargroup,isize,irank, auxvectord, auxvectori)
13489 IF (isize == irank)
THEN
13493 print *,
' CKPGRP ', ipgrp,
globalparlabelindex(1,itgbi), isize, irank,
' rank deficit !!!'
13511 INTEGER(mpl) :: nan
13512 INTEGER(mpl) :: neg
13514 print *,
' Checking global matrix(D) for NANs ',
size(
globalmatd,kind=mpl)
13519 print *,
' i, nan ', i, nan
13525 print *,
' Checking diagonal elements ',
nagb
13530 print *,
' i, neg ', i, neg
13534 print *,
' CHKMAT summary ', nan, neg
13556 REAL(mpd),
INTENT(IN) :: chi2
13557 INTEGER(mpi),
INTENT(IN) :: ithrd
13558 INTEGER(mpi),
INTENT(IN) :: ndf
13559 REAL(mpd),
INTENT(IN) :: dw
13561 INTEGER(mpl) ::nadd
13589 REAL(mpd),
INTENT(OUT) ::chi2
13590 INTEGER(mpl),
INTENT(OUT) ::ndf
13591 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