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