Millepede-II  V04-09-02
pede.f90
Go to the documentation of this file.
1 ! Code converted using TO_F90 by Alan Miller
2 ! Date: 2012-03-16 Time: 11:06:00
3 
27 
198 
362 
592 
627 
629 PROGRAM mptwo
630  USE mpmod
631  USE mpdalc
632  USE mptest1, ONLY: nplan,del,dvd
633  USE mptest2, ONLY: nlyr,nmx,nmy,sdevx,sdevy
634 
635  IMPLICIT NONE
636  REAL(mps) :: andf
637  REAL(mps) :: c2ndf
638  REAL(mps) :: deltat
639  REAL(mps) :: diff
640  REAL(mps) :: err
641  REAL(mps) :: gbu
642  REAL(mps) :: gmati
643  REAL(mps) :: rej
644  REAL :: rloop1
645  REAL :: rloop2
646  REAL :: rstext
647  REAL(mps) :: secnd
648  REAL :: rst
649  REAL :: rstp
650  REAL, DIMENSION(2) :: ta
651  INTEGER(mpi) :: i
652  INTEGER(mpi) :: ii
653  INTEGER(mpi) :: iopnmp
654  INTEGER(mpi) :: ix
655  INTEGER(mpi) :: ixv
656  INTEGER(mpi) :: iy
657  INTEGER(mpi) :: k
658  INTEGER(mpi) :: kfl
659  INTEGER(mpi) :: lun
660  INTEGER :: minut
661  INTEGER :: nhour
662  INTEGER(mpi) :: nmxy
663  INTEGER(mpi) :: nrc
664  INTEGER(mpi) :: nsecnd
665  INTEGER(mpi) :: ntot
666  INTEGER(mpi) :: ntsec
667 
668  CHARACTER (LEN=24) :: chdate
669  CHARACTER (LEN=24) :: chost
670 #ifdef LAPACK64
671  CHARACTER (LEN=6) :: c6
672 #endif
673 
674  INTEGER(mpl) :: rows
675  INTEGER(mpl) :: cols
676 
677  REAL(mpd) :: sums(9)
678  !$ INTEGER(mpi) :: OMP_GET_NUM_PROCS,OMP_GET_MAX_THREADS
679  !$ INTEGER(mpi) :: MXTHRD
680  !$ INTEGER(mpi) :: NPROC
681 
682  REAL etime
683 
684  SAVE
685  ! ...
686  rstp=etime(ta)
687  CALL fdate(chdate)
688 
689  ! millepede monitoring file
690  lunmon=0
691  ! millepede.log file
692  lunlog=8
693  lvllog=1
694  CALL mvopen(lunlog,'millepede.log')
695  CALL getenv('HOSTNAME',chost)
696  IF (chost(1:1) == ' ') CALL getenv('HOST',chost)
697  WRITE(*,*) '($Rev: 210 $)'
698  iopnmp=0
699  !$ iopnmp=1
700  !$ WRITE(*,*) 'using OpenMP (TM)'
701 #ifdef LAPACK64
702  WRITE(*,*) 'using LAPACK64 with ', lapack64
703 #endif
704 #ifdef __GFORTRAN__
705  WRITE(*,111) __gnuc__ , __gnuc_minor__ , __gnuc_patchlevel__
706 111 FORMAT(' compiled with gcc ',i0,'.',i0,'.',i0)
707 #endif
708 #ifdef __PGIC__
709  WRITE(*,111) __pgic__ , __pgic_minor__ , __pgic_patchlevel__
710 111 FORMAT(' compiled with pgi ',i0,'.',i0,'.',i0)
711 #endif
712  WRITE(*,*) ' '
713  WRITE(*,*) ' < Millepede II-P starting ... ',chdate
714  WRITE(*,*) ' ',chost
715  WRITE(*,*) ' '
716 
717  WRITE(8,*) '($Rev: 210 $)'
718  WRITE(8,*) ' '
719  WRITE(8,*) 'Log-file Millepede II-P ', chdate
720  WRITE(8,*) ' ', chost
721 
722  CALL peend(-1,'Still running or crashed')
723  ! read command line and text files
724 
725  CALL filetc ! command line and steering file analysis
726  CALL filetx ! read text files
727 
728  IF (icheck > 0) THEN
729  WRITE(*,*) '!!! Checking input only, no calculation of a solution !!!'
730  WRITE(8,*) '!!! Checking input only, no calculation of a solution !!!'
731  END IF
732  lvllog=mprint ! export print level
733  IF (memdbg > 0) printflagalloc=1 ! debug memory management
734  !$ WRITE(*,*)
735  !$ NPROC=1
736  !$ MXTHRD=1
737  !$ NPROC=OMP_GET_NUM_PROCS() ! number of processors available
738  !$ CALL OMP_SET_NUM_THREADS(MTHRD) ! set max number of threads to MTHRD
739  !$ MXTHRD=OMP_GET_MAX_THREADS() ! get max number of threads back
740  !$ WRITE(*,*) 'Number of processors available: ', NPROC
741  !$ WRITE(*,*) 'Maximum number of OpenMP threads: ', MXTHRD
742  !$ WRITE(*,*) 'Number of threads for processing: ', MTHRD
743  !$ IF (MXREC.GT.0) MTHRDR=1 ! to get allways the same MXREC records
744  !$ IF (ICHECK.GT.1) MTHRDR=1 ! to get allways the same order of records
745  !$ WRITE(*,*) 'Number of threads for reading: ', MTHRDR
746  !$POMP INST INIT ! start profiling with ompP
747 #ifdef LAPACK64
748  IF(iopnmp > 0) THEN
749  CALL getenv('OMP_NUM_THREADS',c6)
750  ELSE
751  CALL getenv(lapack64//'_NUM_THREADS',c6)
752  END IF
753  IF (c6(1:1) == ' ') THEN
754  IF(iopnmp > 0) THEN
755  WRITE(*,*) 'Number of threads for LAPACK: unkown (empty OMP_NUM_THREADS)'
756  ELSE
757  WRITE(*,*) 'Number of threads for LAPACK: unkown (empty ',lapack64//'_NUM_THREADS)'
758  END IF
759  ELSE
760  WRITE(*,*) 'Number of threads for LAPACK: ', c6
761  END IF
762 #endif
763  IF (ncache < 0) THEN
764  ncache=25000000*mthrd ! default cache size (100 MB per thread)
765  ENDIF
766  rows=6; cols=mthrdr
767  CALL mpalloc(readbufferinfo,rows,cols,'read buffer header')
768  ! histogram file
769  lun=7
770  CALL mvopen(lun,'millepede.his')
771  CALL hmplun(lun) ! unit for histograms
772  CALL gmplun(lun) ! unit for xy data
773 
774  ! debugging
775  IF(nrecpr /= 0.OR.nrecp2 /= 0) THEN
776  CALL mvopen(1,'mpdebug.txt')
777  END IF
778 
779  rstext=etime(ta)
780  times(0)=rstext-rstp ! time for text processing
781 
782  ! preparation of data sub-arrays
783 
784  CALL loop1
785  rloop1=etime(ta)
786  times(1)=rloop1-rstext ! time for LOOP1
787 
788  CALL loop2
789  IF(chicut /= 0.0) THEN
790  WRITE(8,*) 'Chi square cut equiv 3 st.dev applied ...'
791  WRITE(8,*) ' in first iteration with factor',chicut
792  WRITE(8,*) ' in second iteration with factor',chirem
793  WRITE(8,*) ' (reduced by sqrt in next iterations)'
794  END IF
795 
796  IF(lhuber /= 0) THEN
797  WRITE(8,*) 'Down-weighting of outliers in', lhuber,' iterations'
798  WRITE(8,*) 'Cut on downweight fraction',dwcut
799  END IF
800 
801  rloop2=etime(ta)
802  times(2)=rloop2-rloop1 ! time for LOOP2
803 
804  IF(icheck > 0) THEN
805  CALL prtstat
806  CALL peend(0,'Ended normally')
807  GOTO 99 ! only checking input
808  END IF
809 
810  ! use different solution methods
811 
812  CALL mstart('Iteration') ! Solution module starting
813 
814  CALL xloopn ! all methods
815 
816  ! ------------------------------------------------------------------
817 
818  IF(nloopn > 2.AND.nhistp /= 0) THEN ! last iteration
819  CALL hmprnt(3) ! scaled residual of single measurement (with global deriv.)
820  CALL hmprnt(12) ! scaled residual of single measurement (no global deriv.)
821  CALL hmprnt(4) ! chi^2/Ndf
822  END IF
823  IF(nloopn > 2) THEN
824  CALL hmpwrt(3)
825  CALL hmpwrt(12)
826  CALL hmpwrt(4)
827  CALL gmpwrt(4) ! location, dispersion (res.) as a function of record nr
828  IF (nloopn <= lfitnp) THEN
829  CALL hmpwrt(13)
830  CALL hmpwrt(14)
831  CALL gmpwrt(5)
832  END IF
833  END IF
834  IF(nhistp /= 0) THEN
835  CALL gmprnt(1)
836  CALL gmprnt(2)
837  END IF
838  CALL gmpwrt(1) ! output of xy data
839  CALL gmpwrt(2) ! output of xy data
840  ! 'track quality' per binary file
841  IF (nfilb > 1) THEN
842  CALL gmpdef(6,1,'log10(#records) vs file number')
843  CALL gmpdef(7,1,'final rejection fraction vs file number')
844  CALL gmpdef(8,1, &
845  'final <Chi^2/Ndf> from accepted local fits vs file number')
846  CALL gmpdef(9,1, '<Ndf> from accepted local fits vs file number')
847 
848  DO i=1,nfilb
849  kfl=kfd(2,i)
850  nrc=-kfd(1,i)
851  IF (nrc > 0) THEN
852  rej=real(nrc-jfd(kfl),mps)/real(nrc,mps)
853  CALL gmpxy(6,real(kfl,mps),log10(real(nrc,mps))) ! log10(#records) vs file
854  CALL gmpxy(7,real(kfl,mps),rej) ! rejection fraction vs file
855  END IF
856  IF (jfd(kfl) > 0) THEN
857  c2ndf=cfd(kfl)/real(jfd(kfl),mps)
858  CALL gmpxy(8,real(kfl,mps),c2ndf) ! <Chi2/NDF> vs file
859  andf=real(dfd(kfl),mps)/real(jfd(kfl),mps)
860  CALL gmpxy(9,real(kfl,mps),andf) ! <NDF> vs file
861  END IF
862  END DO
863  IF(nhistp /= 0) THEN
864  CALL gmprnt(6)
865  CALL gmprnt(7)
866  CALL gmprnt(8)
867  CALL gmprnt(9)
868  END IF
869  CALL gmpwrt(6) ! output of xy data
870  CALL gmpwrt(7) ! output of xy data
871  CALL gmpwrt(8) ! output of xy data
872  CALL gmpwrt(9) ! output of xy data
873  END IF
874 
875  IF(ictest == 1) THEN
876  WRITE(*,*) ' '
877  WRITE(*,*) 'Misalignment test wire chamber'
878  WRITE(*,*) ' '
879 
880  CALL hmpdef( 9,-0.0015,+0.0015,'True - fitted displacement')
881  CALL hmpdef(10,-0.0015,+0.0015,'True - fitted Vdrift')
882  DO i=1,4
883  sums(i)=0.0_mpd
884  END DO
885  DO i=1,nplan
886  diff=real(-del(i)-globalparameter(i),mps)
887  sums(1)=sums(1)+diff
888  sums(2)=sums(2)+diff*diff
889  diff=real(-dvd(i)-globalparameter(100+i),mps)
890  sums(3)=sums(3)+diff
891  sums(4)=sums(4)+diff*diff
892  END DO
893  sums(1)=0.01_mpd*sums(1)
894  sums(2)=sqrt(0.01_mpd*sums(2))
895  sums(3)=0.01_mpd*sums(3)
896  sums(4)=sqrt(0.01_mpd*sums(4))
897  WRITE(*,143) 'Parameters 1 - 100: mean =',sums(1), 'rms =',sums(2)
898  WRITE(*,143) 'Parameters 101 - 200: mean =',sums(3), 'rms =',sums(4)
899 143 FORMAT(6x,a28,f9.6,3x,a5,f9.6)
900  WRITE(*,*) ' '
901  WRITE(*,*) ' '
902  WRITE(*,*) ' I '
903  WRITE(*,*) ' --- '
904  DO i=1,100
905  WRITE(*,102) i,-del(i),globalparameter(i),-del(i)-globalparameter(i), &
906  -dvd(i),globalparameter(100+i),-dvd(i)-globalparameter(100+i)
907  diff=real(-del(i)-globalparameter(i),mps)
908  CALL hmpent( 9,diff)
909  diff=real(-dvd(i)-globalparameter(100+i),mps)
910  CALL hmpent(10,diff)
911  END DO
912  IF(nhistp /= 0) THEN
913  CALL hmprnt( 9)
914  CALL hmprnt(10)
915  END IF
916  CALL hmpwrt( 9)
917  CALL hmpwrt(10)
918  END IF
919  IF(ictest > 1) THEN
920  WRITE(*,*) ' '
921  WRITE(*,*) 'Misalignment test Si tracker'
922  WRITE(*,*) ' '
923 
924  CALL hmpdef( 9,-0.0025,+0.0025,'True - fitted displacement X')
925  CALL hmpdef(10,-0.025,+0.025,'True - fitted displacement Y')
926  DO i=1,9
927  sums(i)=0.0_mpd
928  END DO
929  nmxy=nmx*nmy
930  ix=0
931  iy=ntot
932  DO i=1,nlyr
933  DO k=1,nmxy
934  ix=ix+1
935  diff=real(-sdevx((i-1)*nmxy+k)-globalparameter(ix),mps)
936  sums(1)=sums(1)+1.0_mpd
937  sums(2)=sums(2)+diff
938  sums(3)=sums(3)+diff*diff
939  ixv=globalparlabelindex(2,ix)
940  IF (ixv > 0.AND.metsol == 1.OR.metsol == 2) THEN
941  ii=(ixv*ixv+ixv)/2
942  gmati=real(globalmatd(ii),mps)
943  err=sqrt(abs(gmati))
944  diff=diff/err
945  sums(7)=sums(7)+1.0_mpd
946  sums(8)=sums(8)+diff
947  sums(9)=sums(9)+diff*diff
948  END IF
949  END DO
950  IF (mod(i,3) == 1) THEN
951  DO k=1,nmxy
952  iy=iy+1
953  diff=-real(sdevy((i-1)*nmxy+k)-globalparameter(iy),mps)
954  sums(4)=sums(4)+1.0_mpd
955  sums(5)=sums(5)+diff
956  sums(6)=sums(6)+diff*diff
957  ixv=globalparlabelindex(2,iy)
958  IF (ixv > 0.AND.metsol == 1.OR.metsol == 2) THEN
959  ii=(ixv*ixv+ixv)/2
960  gmati=real(globalmatd(ii),mps)
961  err=sqrt(abs(gmati))
962  diff=diff/err
963  sums(7)=sums(7)+1.0_mpd
964  sums(8)=sums(8)+diff
965  sums(9)=sums(9)+diff*diff
966  END IF
967  END DO
968  END IF
969  END DO
970  sums(2)=sums(2)/sums(1)
971  sums(3)=sqrt(sums(3)/sums(1))
972  sums(5)=sums(5)/sums(4)
973  sums(6)=sqrt(sums(6)/sums(4))
974  WRITE(*,143) 'Parameters 1 - 500: mean =',sums(2), 'rms =',sums(3)
975  WRITE(*,143) 'Parameters 501 - 700: mean =',sums(5), 'rms =',sums(6)
976  IF (sums(7) > 0.5_mpd) THEN
977  sums(8)=sums(8)/sums(7)
978  sums(9)=sqrt(sums(9)/sums(7))
979  WRITE(*,143) 'Parameter pulls, all: mean =',sums(8), 'rms =',sums(9)
980  END IF
981  WRITE(*,*) ' '
982  WRITE(*,*) ' '
983  WRITE(*,*) ' I '
984  WRITE(*,*) ' --- '
985  ix=0
986  iy=ntot
987  DO i=1,nlyr
988  DO k=1,nmxy
989  ix=ix+1
990  diff=real(-sdevx((i-1)*nmxy+k)-globalparameter(ix),mps)
991  CALL hmpent( 9,diff)
992  WRITE(*,102) ix,-sdevx((i-1)*nmxy+k),globalparameter(ix),-diff
993  END DO
994  END DO
995  DO i=1,nlyr
996  IF (mod(i,3) == 1) THEN
997  DO k=1,nmxy
998  iy=iy+1
999  diff=real(-sdevy((i-1)*nmxy+k)-globalparameter(iy),mps)
1000  CALL hmpent(10,diff)
1001  WRITE(*,102) iy,-sdevy((i-1)*nmxy+k),globalparameter(iy),-diff
1002  END DO
1003  END IF
1004  END DO
1005  IF(nhistp /= 0) THEN
1006  CALL hmprnt( 9)
1007  CALL hmprnt(10)
1008  END IF
1009  CALL hmpwrt( 9)
1010  CALL hmpwrt(10)
1011  END IF
1012 
1013  IF(nrec1+nrec2 > 0) THEN
1014  WRITE(8,*) ' '
1015  IF(nrec1 > 0) THEN
1016  WRITE(8,*) 'Record',nrec1,' has largest residual:',value1
1017  END IF
1018  IF(nrec2 > 0) THEN
1019  WRITE(8,*) 'Record',nrec2,' has largest Chi^2/Ndf:',value2
1020  END IF
1021  END IF
1022  IF(nrec3 < huge(nrec3)) THEN
1023  WRITE(8,*) 'Record',nrec3, ' is first with error (rank deficit/NaN)'
1024  END IF
1025 99 WRITE(8,*) ' '
1026  IF (iteren > mreqenf) THEN
1027  WRITE(8,*) 'In total 3 +',nloopn,' loops through the data files'
1028  ELSE
1029  WRITE(8,*) 'In total 2 +',nloopn,' loops through the data files'
1030  ENDIF
1031  IF (mnrsit > 0) THEN
1032  WRITE(8,*) ' '
1033  WRITE(8,*) 'In total ',mnrsit,' internal MINRES iterations'
1034  END IF
1035 
1036  WRITE(8,103) times(0),times(1),times(2),times(4),times(7), &
1037  times(5),times(8),times(3),times(6)
1038 
1039  rst=etime(ta)
1040  deltat=rst-rstp
1041  ntsec=nint(deltat,mpi)
1042  CALL sechms(deltat,nhour,minut,secnd)
1043  nsecnd=nint(secnd,mpi) ! round
1044  WRITE(8,*) 'Total time =',ntsec,' seconds =',nhour,' h',minut, &
1045  ' m',nsecnd,' seconds'
1046  CALL fdate(chdate)
1047  WRITE(8,*) 'end ', chdate
1048  gbu=1.0e-9*real(maxwordsalloc*(bit_size(1_mpi)/8),mps) ! GB used
1049  WRITE(8,*) ' '
1050  WRITE(8,105) gbu
1051 
1052  ! Rejects ----------------------------------------------------------
1053 
1054  IF(nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) /= 0) THEN
1055  WRITE(8,*) ' '
1056  WRITE(8,*) 'Data rejected in last iteration: '
1057  WRITE(8,*) ' ', &
1058  nrejec(0), ' (rank deficit/NaN) ',nrejec(1),' (Ndf=0) ', &
1059  nrejec(2), ' (huge) ',nrejec(3),' (large)'
1060  WRITE(8,*) ' '
1061  END IF
1062  IF (icheck <= 0) CALL explfc(8)
1063 
1064  WRITE(*,*) ' '
1065  WRITE(*,*) ' < Millepede II-P ending ... ', chdate ! with exit code',ITEXIT,' >'
1066  WRITE(*,*) ' '
1067  gbu=1.0e-9*real(maxwordsalloc*(bit_size(1_mpi)/8),mps) ! GB used
1068  WRITE(*,105) gbu
1069  WRITE(*,*) ' '
1070 
1071 102 FORMAT(2x,i4,2x,3f10.5,2x,3f10.5)
1072 103 FORMAT(' Times [in sec] for text processing',f12.3/ &
1073  ' LOOP1',f12.3/ &
1074  ' LOOP2',f12.3/ &
1075  ' func. value ',f12.3,' *',f4.0/ &
1076  ' func. value, global matrix, solution',f12.3,' *',f4.0/ &
1077  ' new solution',f12.3,' *',f4.0/)
1078 105 FORMAT(' Peak dynamic memory allocation: ',f11.6,' GB')
1079 END PROGRAM mptwo ! Mille
1080 
1087 
1088 SUBROUTINE solglo(ivgbi)
1089  USE mpmod
1090  USE minresmodule, ONLY: minres
1091 
1092  IMPLICIT NONE
1093  REAL(mps) :: par
1094  REAL(mps) :: dpa
1095  REAL(mps) :: err
1096  REAL(mps) :: gcor2
1097  INTEGER(mpi) :: iph
1098  INTEGER(mpi) :: istop
1099  INTEGER(mpi) :: itgbi
1100  INTEGER(mpi) :: itgbl
1101  INTEGER(mpi) :: itn
1102  INTEGER(mpi) :: itnlim
1103  INTEGER(mpi) :: nout
1104 
1105  INTEGER(mpi), INTENT(IN) :: ivgbi
1106 
1107  REAL(mpd) :: shift
1108  REAL(mpd) :: rtol
1109  REAL(mpd) :: anorm
1110  REAL(mpd) :: acond
1111  REAL(mpd) :: arnorm
1112  REAL(mpd) :: rnorm
1113  REAL(mpd) :: ynorm
1114  REAL(mpd) :: gmati
1115  REAL(mpd) :: diag
1116  REAL(mpd) :: matij
1117  LOGICAL :: checka
1118  EXTERNAL avprod, mcsolv, mvsolv
1119  SAVE
1120  DATA iph/0/
1121  ! ...
1122  IF(iph == 0) THEN
1123  iph=1
1124  WRITE(*,101)
1125  END IF
1126  itgbi=globalparvartototal(ivgbi)
1127  itgbl=globalparlabelindex(1,itgbi)
1128 
1129  globalvector=0.0_mpd ! reset rhs vector IGVEC
1130  globalvector(ivgbi)=1.0_mpd
1131 
1132  ! NOUT =6
1133  nout =0
1134  itnlim=200
1135  shift =0.0_mpd
1136  rtol = mrestl ! from steering
1137  checka=.false.
1138 
1139 
1140  IF(mbandw == 0) THEN ! default preconditioner
1141  CALL minres(nagb, avprod, mcsolv, globalvector, shift, checka ,.true. , &
1142  globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1143 
1144  ELSE IF(mbandw > 0) THEN ! band matrix preconditioner
1145  CALL minres(nagb, avprod, mvsolv, globalvector, shift, checka ,.true. , &
1146  globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1147  ELSE
1148  CALL minres(nagb, avprod, mvsolv, globalvector, shift, checka ,.false. , &
1149  globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1150  END IF
1151 
1152  par=real(globalparameter(itgbi),mps)
1153  dpa=real(par-globalparstart(itgbi),mps)
1154  gmati=globalcorrections(ivgbi)
1155  err=sqrt(abs(real(gmati,mps)))
1156  IF(gmati < 0.0_mpd) err=-err
1157  diag=matij(ivgbi,ivgbi)
1158  gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps) ! global correlation (squared)
1159  WRITE(*,102) itgbl,par,real(globalparpresigma(itgbi),mps),dpa,err,gcor2,itn
1160 101 FORMAT(1x,' label parameter presigma differ', &
1161  ' Error gcor^2 iit'/ 1x,'---------',2x,5('-----------'),2x,'----')
1162 102 FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1163 END SUBROUTINE solglo
1164 
1171 
1172 SUBROUTINE solgloqlp(ivgbi)
1173  USE mpmod
1174  USE minresqlpmodule, ONLY: minresqlp
1175 
1176  IMPLICIT NONE
1177  REAL(mps) :: par
1178  REAL(mps) :: dpa
1179  REAL(mps) :: err
1180  REAL(mps) :: gcor2
1181  INTEGER(mpi) :: iph
1182  INTEGER(mpi) :: istop
1183  INTEGER(mpi) :: itgbi
1184  INTEGER(mpi) :: itgbl
1185  INTEGER(mpi) :: itn
1186  INTEGER(mpi) :: itnlim
1187  INTEGER(mpi) :: nout
1188 
1189  INTEGER(mpi), INTENT(IN) :: ivgbi
1190 
1191  REAL(mpd) :: shift
1192  REAL(mpd) :: rtol
1193  REAL(mpd) :: mxxnrm
1194  REAL(mpd) :: trcond
1195  REAL(mpd) :: gmati
1196  REAL(mpd) :: diag
1197  REAL(mpd) :: matij
1198 
1199  EXTERNAL avprod, mcsolv, mvsolv
1200  SAVE
1201  DATA iph/0/
1202  ! ...
1203  IF(iph == 0) THEN
1204  iph=1
1205  WRITE(*,101)
1206  END IF
1207  itgbi=globalparvartototal(ivgbi)
1208  itgbl=globalparlabelindex(1,itgbi)
1209 
1210  globalvector=0.0_mpd ! reset rhs vector IGVEC
1211  globalvector(ivgbi)=1.0_mpd
1212 
1213  ! NOUT =6
1214  nout =0
1215  itnlim=200
1216  shift =0.0_mpd
1217  rtol = mrestl ! from steering
1218  mxxnrm = real(nagb,mpd)/sqrt(epsilon(mxxnrm))
1219  IF(mrmode == 1) THEN
1220  trcond = 1.0_mpd/epsilon(trcond) ! only QR
1221  ELSE IF(mrmode == 2) THEN
1222  trcond = 1.0_mpd ! only QLP
1223  ELSE
1224  trcond = mrtcnd ! QR followed by QLP
1225  END IF
1226 
1227  IF(mbandw == 0) THEN ! default preconditioner
1228  CALL minresqlp( n=nagb, aprod=avprod, b=globalvector, msolve=mcsolv, nout=nout, &
1229  itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1230  x=globalcorrections, istop=istop, itn=itn)
1231  ELSE IF(mbandw > 0) THEN ! band matrix preconditioner
1232  CALL minresqlp( n=nagb, aprod=avprod, b=globalvector, msolve=mvsolv, nout=nout, &
1233  itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1234  x=globalcorrections, istop=istop, itn=itn)
1235  ELSE
1236  CALL minresqlp( n=nagb, aprod=avprod, b=globalvector, nout=nout, &
1237  itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1238  x=globalcorrections, istop=istop, itn=itn)
1239  END IF
1240 
1241  par=real(globalparameter(itgbi),mps)
1242  dpa=real(par-globalparstart(itgbi),mps)
1243  gmati=globalcorrections(ivgbi)
1244  err=sqrt(abs(real(gmati,mps)))
1245  IF(gmati < 0.0_mpd) err=-err
1246  diag=matij(ivgbi,ivgbi)
1247  gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps) ! global correlation (squared)
1248  WRITE(*,102) itgbl,par,real(globalparpresigma(itgbi),mps),dpa,err,gcor2,itn
1249 101 FORMAT(1x,' label parameter presigma differ', &
1250  ' Error gcor^2 iit'/ 1x,'---------',2x,5('-----------'),2x,'----')
1251 102 FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1252 END SUBROUTINE solgloqlp
1253 
1255 SUBROUTINE addcst
1256  USE mpmod
1257 
1258  IMPLICIT NONE
1259  REAL(mpd) :: climit
1260  REAL(mpd) :: factr
1261  REAL(mpd) :: sgm
1262 
1263  INTEGER(mpi) :: i
1264  INTEGER(mpi) :: icgb
1265  INTEGER(mpi) :: irhs
1266  INTEGER(mpi) :: itgbi
1267  INTEGER(mpi) :: ivgb
1268  INTEGER(mpi) :: j
1269  INTEGER(mpi) :: jcgb
1270  INTEGER(mpi) :: l
1271  INTEGER(mpi) :: label
1272  INTEGER(mpi) :: nop
1273  INTEGER(mpi) :: inone
1274 
1275  REAL(mpd) :: rhs
1276  REAL(mpd) :: drhs(4)
1277  INTEGER(mpi) :: idrh (4)
1278  SAVE
1279  ! ...
1280  nop=0
1281  IF(lenconstraints == 0) RETURN ! no constraints
1282  climit=1.0e-5 ! limit for printout
1283  irhs=0 ! number of values in DRHS(.), to be printed
1284 
1285  DO jcgb=1,ncgb
1286  icgb=matconssort(3,jcgb) ! unsorted constraint index
1287  i=vecconsstart(icgb)
1288  rhs=listconstraints(i )%value ! right hand side
1289  sgm=listconstraints(i+1)%value ! sigma parameter
1290  DO j=i+2,vecconsstart(icgb+1)-1
1291  label=listconstraints(j)%label
1292  factr=listconstraints(j)%value
1293  itgbi=inone(label) ! -> ITGBI= index of parameter label
1294  ivgb =globalparlabelindex(2,itgbi) ! -> variable-parameter index
1295 
1296  IF(icalcm == 1.AND.nagb > nvgb.AND.ivgb > 0) THEN
1297  CALL mupdat(nvgb+jcgb,ivgb,factr) ! add to matrix
1298  END IF
1299 
1300  rhs=rhs-factr*globalparameter(itgbi) ! reduce residuum
1301  END DO
1302  IF(abs(rhs) > climit) THEN
1303  irhs=irhs+1
1304  idrh(irhs)=jcgb
1305  drhs(irhs)=rhs
1306  nop=1
1307  IF(irhs == 4) THEN
1308  WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1309  irhs=0
1310  END IF
1311  END IF
1312  vecconsresiduals(jcgb)=rhs
1313  IF (nagb > nvgb) globalvector(nvgb+jcgb)=rhs
1314  END DO
1315 
1316  IF(irhs /= 0) THEN
1317  WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1318  END IF
1319  IF(nop == 0) RETURN
1320  WRITE(*,102) ' Constraints: only equation values >', climit,' are printed'
1321 101 FORMAT(' ',4(i4,g11.3))
1322 102 FORMAT(a,g11.2,a)
1323 END SUBROUTINE addcst
1324 
1328 
1329 SUBROUTINE prpcon
1330  USE mpmod
1331  USE mpdalc
1332 
1333  IMPLICIT NONE
1334  INTEGER(mpi) :: i
1335  INTEGER(mpi) :: icgb
1336  INTEGER(mpi) :: isblck
1337  INTEGER(mpi) :: ifrst
1338  INTEGER(mpi) :: ilast
1339  INTEGER(mpi) :: itgbi
1340  INTEGER(mpi) :: ivgb
1341  INTEGER(mpi) :: jcgb
1342  INTEGER(mpi) :: label
1343  INTEGER(mpi) :: labelf
1344  INTEGER(mpi) :: labell
1345  INTEGER(mpi) :: ncon
1346  INTEGER(mpi) :: npar
1347  INTEGER(mpi) :: nconmx
1348  INTEGER(mpi) :: nparmx
1349  INTEGER(mpi) :: inone
1350  INTEGER(mpi) :: itype
1351  INTEGER(mpi) :: ncgbw
1352  INTEGER(mpi) :: newlen
1353  INTEGER(mpi) :: nvar
1354  INTEGER(mpi) :: last
1355  INTEGER(mpi) :: lastlen
1356 
1357  INTEGER(mpl):: length
1358  INTEGER(mpl) :: rows
1359 
1360  ncgb=0
1361  ncgbw=0
1362  ncgbe=0
1363  IF(lenconstraints == 0) RETURN ! no constraints
1364 
1365  newlen=0
1366  lastlen=0
1367  nvar=-1
1368  i=0
1369  last=-1
1370  itype=0
1371  ! find next constraint header and count nr of constraints
1372  DO WHILE(i < lenconstraints)
1373  i=i+1
1374  label=listconstraints(i)%label
1375  IF(last == 0.AND.label < 0) THEN
1376  IF (ncgb > 0 .AND. icheck>0) WRITE(*,113) ncgb, newlen-lastlen-3, nvar
1377  IF (nvar == 0) ncgbe=ncgbe+1
1378  IF (nvar == 0 .AND. iskpec > 0) THEN
1379  ! overwrite
1380  newlen=lastlen
1381  ! copy previous value (for label 0)
1382  newlen=newlen+1
1383  listconstraints(newlen)%value=listconstraints(i-1)%value
1384  ELSE
1385  lastlen=newlen-1 ! end of last accepted constraint
1386  END IF
1387  ncgb=ncgb+1
1388  itype=-label
1389  IF(itype == 2) ncgbw=ncgbw+1
1390  nvar=0
1391  END IF
1392  last=label
1393  IF(label > 0) THEN
1394  itgbi=inone(label) ! -> ITGBI= index of parameter label
1395  ivgb =globalparlabelindex(2,itgbi) ! -> variable-parameter index
1396  IF (ivgb > 0) nvar=nvar+1
1397  END IF
1398  IF(label > 0.AND.itype == 2) THEN ! weighted constraints
1399  itgbi=inone(label) ! -> ITGBI= index of parameter label
1400  listconstraints(i)%value=listconstraints(i)%value*globalparcounts(itgbi)
1401  END IF
1402  newlen=newlen+1
1403  listconstraints(newlen)%label=listconstraints(i)%label ! copy label
1404  listconstraints(newlen)%value=listconstraints(i)%value ! copy value
1405  END DO
1406  IF (ncgb > 0 .AND. icheck>0) WRITE(*,113) ncgb, newlen-lastlen-2, nvar
1407  IF (nvar == 0) ncgbe=ncgbe+1
1408  IF (nvar == 0 .AND. iskpec > 0) newlen=lastlen
1409  lenconstraints=newlen
1410 
1411  IF (ncgbe > 0 .AND. iskpec > 0) THEN
1412  WRITE(*,*) 'PRPCON:',ncgbe,' empty constraints skipped'
1413  ncgb=ncgb-ncgbe
1414  END IF
1415  IF (ncgbw == 0) THEN
1416  WRITE(*,*) 'PRPCON:',ncgb,' constraints accepted'
1417  ELSE
1418  WRITE(*,*) 'PRPCON:',ncgb,' constraints accepted,',ncgbw, 'weighted'
1419  END IF
1420  WRITE(*,*)
1421 
1422  IF(lenconstraints == 0) RETURN ! no constraints left
1423 
1424  ! keys and index for sorting of constraints
1425  length=ncgb+1; rows=3
1426  CALL mpalloc(matconssort,rows,length,'keys and index for sorting (I)')
1427  matconssort(1,ncgb+1)=ntgb+1
1428  ! start of constraint in list
1429  CALL mpalloc(vecconsstart,length,'start of constraint in list (I)')
1431  ! start and parameter range of constraint blocks
1432  CALL mpalloc(matconsblocks,rows,length,'start of constraint blocks, par. range (I)')
1433 
1434  ! prepare
1435  i=1
1436  DO icgb=1,ncgb
1437  ! new constraint
1438  vecconsstart(icgb)=i
1439  matconssort(1,icgb)=ntgb ! min variable parameter
1440  matconssort(2,icgb)=0 ! max variable parameter
1441  matconssort(3,icgb)=icgb ! index
1442  i=i+2
1443  DO
1444  label=listconstraints(i)%label
1445  itgbi=inone(label) ! -> ITGBI= index of parameter label
1446  ivgb =globalparlabelindex(2,itgbi) ! -> variable-parameter index
1447  IF(ivgb > 0) THEN
1448  matconssort(1,icgb)=min(matconssort(1,icgb),ivgb)
1449  matconssort(2,icgb)=max(matconssort(2,icgb),ivgb)
1450  END IF
1451  i=i+1
1452  IF(i > lenconstraints) EXIT
1453  IF(listconstraints(i)%label == 0) EXIT
1454  END DO
1455  END DO
1456  ! force single 'full size' block
1457  !matConsSort(1,1)=1
1458  !matConsSort(2,1)=nvgb
1459  ! sort constraints
1460  CALL sort2i(matconssort,ncgb)
1461 
1462  ! loop over sorted constraints, try to split into blocks
1463  ncblck=0
1464  nconmx=0
1465  nparmx=0
1466  mszcon=0
1467  mszprd=0
1468  isblck=1
1469  ilast=0
1470  DO jcgb=1,ncgb
1471  ! index in list
1472  icgb=matconssort(3,jcgb)
1473  ! split into disjoint blocks
1474  ilast=max(ilast, matconssort(2,jcgb))
1475  IF (icheck > 1) THEN
1476  IF (matconssort(2,jcgb) > matconssort(1,jcgb)) THEN
1479  ELSE
1480  labelf=0; labell=0
1481  END IF
1482  WRITE(*,*) ' Cons. sorted', jcgb, icgb, vecconsstart(icgb), labelf, labell
1483  END IF
1484  IF (matconssort(1,jcgb+1) > ilast) THEN
1485  ncblck=ncblck+1
1486  ifrst=matconssort(1,isblck)
1487  IF (ifrst > ilast) ifrst=ilast+1 ! empty constraint block (enforce npar=0)
1488  matconsblocks(1,ncblck)=isblck
1489  matconsblocks(2,ncblck)=ifrst ! save first parameter in block
1490  matconsblocks(3,ncblck)=ilast ! save last parameter in block
1491  ncon=jcgb+1-isblck
1492  npar=ilast+1-ifrst
1493  nconmx=max(nconmx,ncon)
1494  nparmx=max(nparmx,npar)
1495  mszcon=mszcon+int(ncon,mpl)*int(npar,mpl) ! (sum of) block size for constraint matrix
1496  mszprd=mszprd+int(ncon,mpl)*int(ncon+1,mpl)/2 ! (sum of) block size for product matrix
1497  IF (icheck > 0) THEN
1498  IF (ilast > ifrst) THEN
1499  labelf=globalparlabelindex(1,globalparvartototal(ifrst))
1500  labell=globalparlabelindex(1,globalparvartototal(ilast))
1501  ELSE
1502  labelf=0; labell=0
1503  END IF
1504  WRITE(*,*) ' Cons. block ', ncblck, isblck, jcgb, labelf, labell
1505  ENDIF
1506  ! reset for new block
1507  isblck=jcgb+1
1508  ! update index ranges
1509  globalindexranges(ifrst)=max(globalindexranges(ifrst),ilast)
1510  END IF
1511  END DO
1512  matconsblocks(1,ncblck+1)=ncgb+1
1513 
1514  IF (ncblck+icheck > 1) THEN
1515  WRITE(*,*)
1516  WRITE(*,*) 'PRPCON: constraints split into ', ncblck, '(disjoint) blocks'
1517  WRITE(*,*) ' max block size (cons., par.) ', nconmx, nparmx
1518  IF (icheck > 0) WRITE(*,*) ' total block matrix sizes ', mszcon, mszprd
1519  END IF
1520 113 FORMAT(' constraint',i6,' : ',i9,' parameters,',i9,' variable')
1521 
1522 END SUBROUTINE prpcon
1523 
1527 
1528 SUBROUTINE feasma
1529  USE mpmod
1530  USE mpdalc
1531 
1532  IMPLICIT NONE
1533  REAL(mpd) :: factr
1534  REAL(mpd) :: sgm
1535  INTEGER(mpi) :: i
1536  INTEGER(mpi) :: iblck
1537  INTEGER(mpi) :: icgb
1538  INTEGER(mpl) :: ij
1539  INTEGER(mpi) :: ifirst
1540  INTEGER(mpi) :: ilast
1541  INTEGER(mpl) :: ioffc
1542  INTEGER(mpl) :: ioffp
1543  INTEGER(mpi) :: irank
1544  INTEGER(mpi) :: ipar0
1545  INTEGER(mpi) :: itgbi
1546  INTEGER(mpi) :: ivgb
1547  INTEGER(mpi) :: j
1548  INTEGER(mpi) :: jcgb
1549  INTEGER(mpl) :: ll
1550  INTEGER(mpi) :: label
1551  INTEGER(mpi) :: ncon
1552  INTEGER(mpi) :: npar
1553  INTEGER(mpi) :: nrank
1554  INTEGER(mpi) :: inone
1555 
1556  REAL(mpd):: rhs
1557  REAL(mpd):: evmax
1558  REAL(mpd):: evmin
1559  INTEGER(mpl):: length
1560  REAL(mpd), DIMENSION(:), ALLOCATABLE :: matConstraintsT
1561  REAL(mpd), DIMENSION(:), ALLOCATABLE :: auxVectorD
1562  INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: auxVectorI
1563  SAVE
1564  ! ...
1565 
1566  IF(ncgb == 0) RETURN ! no constraints
1567 
1568  ! product matrix A A^T (A is stored as transposed)
1569  length=mszprd
1570  CALL mpalloc(matconsproduct, length, 'product matrix of constraints (blocks)')
1571  matconsproduct=0.0_mpd
1572  length=ncgb
1573  CALL mpalloc(vecconsresiduals, length, 'residuals of constraints')
1574  CALL mpalloc(vecconssolution, length, 'solution for constraints')
1575  CALL mpalloc(auxvectori,length,'auxiliary array (I)') ! int aux 1
1576  CALL mpalloc(auxvectord,length,'auxiliary array (D)') ! double aux 1
1577  ! constraint matrix A (A is stored as transposed)
1578  length = mszcon
1579  CALL mpalloc(matconstraintst,length,'transposed matrix of constraints (blocks)')
1580  matconstraintst=0.0_mpd
1581 
1582  ! loop over sorted constraints, fill matrices, get rank, inverted product matrix (in blocks)
1583  ioffc=0 ! block offset in constraint matrix
1584  ioffp=0 ! block offset in product matrix
1585  nrank=0
1586  DO iblck=1,ncblck
1587  ifirst=matconsblocks(1,iblck) ! first constraint in block
1588  ilast=matconsblocks(1,iblck+1)-1 ! last constraint in block
1589  ncon=ilast+1-ifirst
1590  ipar0=matconsblocks(2,iblck)-1 ! parameter offset
1591  npar=matconsblocks(3,iblck)-ipar0 ! number of parameters
1592  DO jcgb=ifirst,ilast
1593  ! index in list
1594  icgb=matconssort(3,jcgb)
1595  ! fill constraint matrix
1596  i=vecconsstart(icgb)
1597  rhs=listconstraints(i )%value ! right hand side
1598  sgm=listconstraints(i+1)%value ! sigma parameter
1599  DO j=i+2,vecconsstart(icgb+1)-1
1600  label=listconstraints(j)%label
1601  factr=listconstraints(j)%value
1602  itgbi=inone(label) ! -> ITGBI= index of parameter label
1603  ivgb =globalparlabelindex(2,itgbi) ! -> variable-parameter index
1604  IF(ivgb > 0) matconstraintst(int(jcgb-ifirst,mpl)*int(npar,mpl)+ivgb-ipar0+ioffc)=factr ! matrix element
1605  globalparcons(itgbi)=globalparcons(itgbi)+1
1606  rhs=rhs-factr*globalparameter(itgbi) ! reduce residuum
1607  END DO
1608  vecconsresiduals(jcgb)=rhs ! constraint discrepancy
1609  END DO
1610 
1611  ! get rank of blocks
1612  DO ll=ioffc+1,ioffc+npar
1613  ij=ioffp
1614  DO i=1,ncon
1615  DO j=1,i
1616  ij=ij+1
1617  matconsproduct(ij)=matconsproduct(ij)+ &
1618  matconstraintst(int(i-1,mpl)*int(npar,mpl)+ll)* &
1619  matconstraintst(int(j-1,mpl)*int(npar,mpl)+ll)
1620  END DO
1621  END DO
1622  END DO
1623  ! inversion of product matrix of constraints
1624  CALL sqminv(matconsproduct(ioffp+1:ij),vecconsresiduals(ifirst:ilast),ncon,irank, auxvectord, auxvectori)
1625  IF (icheck > 1) WRITE(*,*) ' Constraint block rank', iblck, ncon, irank
1626  nrank=nrank+irank
1627  ioffc=ioffc+int(npar,mpl)*int(ncon,mpl)
1628  ioffp=ij
1629  END DO
1630 
1631  nmiss1=ncgb-nrank
1632 
1633  WRITE(*,*) ' '
1634  WRITE(*,*) 'Rank of product matrix of constraints is',nrank, &
1635  ' for',ncgb,' constraint equations'
1636  WRITE(8,*) 'Rank of product matrix of constraints is',nrank, &
1637  ' for',ncgb,' constraint equations'
1638  IF(nrank < ncgb) THEN
1639  WRITE(*,*) 'Warning: insufficient constraint equations!'
1640  WRITE(8,*) 'Warning: insufficient constraint equations!'
1641  IF (iforce == 0) THEN
1642  isubit=1
1643  WRITE(*,*) ' --> enforcing SUBITO mode'
1644  WRITE(8,*) ' --> enforcing SUBITO mode'
1645  END IF
1646  END IF
1647 
1648  ! QL decomposition
1649  IF (nfgb < nvgb) THEN
1650  print *
1651  print *, 'QL decomposition of constraints matrix'
1652  ! monitor progress
1653  IF(monpg1 > 0) THEN
1654  WRITE(lunlog,*) 'QL decomposition of constraints matrix'
1655  CALL monini(lunlog,monpg1,monpg2)
1656  END IF
1657  IF(matsto > 0) THEN ! True unless LAPACK
1658  ! QL decomposition
1659  CALL qlini(nvgb,ncgb,npblck,monpg1)
1660  ! loop over parameter blocks
1661  CALL qldecb(matconstraintst,matparblockoffsets,matconsblocks)
1662  ! check eignevalues of L
1663  CALL qlgete(evmin,evmax)
1664 #ifdef LAPACK64
1665  ELSE
1666  CALL lpqldec(matconstraintst,evmin,evmax)
1667 #endif
1668  END IF
1669  IF(monpg1 > 0) CALL monend()
1670  print *, ' largest |eigenvalue| of L: ', evmax
1671  print *, ' smallest |eigenvalue| of L: ', evmin
1672  IF (evmin == 0.0_mpd) THEN
1673  CALL peend(27,'Aborted, singular QL decomposition of constraints matrix')
1674  stop 'FEASMA: stopping due to singular QL decomposition of constraints matrix'
1675  END IF
1676  END IF
1677 
1678  CALL mpdealloc(matconstraintst)
1679  CALL mpdealloc(auxvectord)
1680  CALL mpdealloc(auxvectori)
1681 
1682  RETURN
1683 END SUBROUTINE feasma ! matrix for feasible solution
1684 
1692 SUBROUTINE feasib(concut,iact)
1693  USE mpmod
1694  USE mpdalc
1695 
1696  IMPLICIT NONE
1697  REAL(mpd) :: factr
1698  REAL(mpd) :: sgm
1699  INTEGER(mpi) :: i
1700  INTEGER(mpi) :: icgb
1701  INTEGER(mpi) :: iter
1702  INTEGER(mpi) :: itgbi
1703  INTEGER(mpi) :: ivgb
1704  INTEGER(mpi) :: iblck
1705  INTEGER(mpi) :: ieblck
1706  INTEGER(mpi) :: isblck
1707  INTEGER(mpi) :: ifirst
1708  INTEGER(mpi) :: ilast
1709  INTEGER(mpi) :: j
1710  INTEGER(mpi) :: jcgb
1711  INTEGER(mpi) :: label
1712  INTEGER(mpi) :: inone
1713  INTEGER(mpi) :: ncon
1714 
1715  REAL(mps), INTENT(IN) :: concut
1716  INTEGER(mpi), INTENT(OUT) :: iact
1717 
1718  REAL(mpd) :: rhs
1719  REAL(mpd) ::sum1
1720  REAL(mpd) ::sum2
1721  REAL(mpd) ::sum3
1722 
1723  REAL(mpd), DIMENSION(:), ALLOCATABLE :: vecCorrections
1724  SAVE
1725 
1726  iact=0
1727  IF(lenconstraints == 0) RETURN ! no constraints
1728 
1729  DO iter=1,2
1730  vecconsresiduals=0.0_mpd
1731 
1732  ! calculate right constraint equation discrepancies
1733  DO jcgb=1,ncgb
1734  icgb=matconssort(3,jcgb) ! unsorted constraint index
1735  i=vecconsstart(icgb)
1736  rhs=listconstraints(i )%value ! right hand side
1737  sgm=listconstraints(i+1)%value ! sigma parameter
1738  DO j=i+2,vecconsstart(icgb+1)-1
1739  label=listconstraints(j)%label
1740  factr=listconstraints(j)%value
1741  itgbi=inone(label) ! -> ITGBI= index of parameter label
1742  rhs=rhs-factr*globalparameter(itgbi) ! reduce residuum
1743  ENDDO
1744  vecconsresiduals(jcgb)=rhs ! constraint discrepancy
1745  END DO
1746 
1747  ! constraint equation discrepancies -------------------------------
1748 
1749  sum1=0.0_mpd
1750  sum2=0.0_mpd
1751  sum3=0.0_mpd
1752  DO icgb=1,ncgb
1753  sum1=sum1+vecconsresiduals(icgb)**2
1754  sum2=sum2+abs(vecconsresiduals(icgb))
1755  sum3=max(sum3,abs(vecconsresiduals(icgb)))
1756  END DO
1757  sum1=sqrt(sum1/real(ncgb,mpd))
1758  sum2=sum2/real(ncgb,mpd)
1759 
1760  IF(iter == 1.AND.sum1 < concut) RETURN ! do nothing if correction small
1761 
1762  IF(iter == 1.AND.ncgb <= 12) THEN
1763  WRITE(*,*) ' '
1764  WRITE(*,*) 'Constraint equation discrepancies:'
1765  WRITE(*,101) (icgb,vecconsresiduals(icgb),icgb=1,ncgb)
1766 101 FORMAT(4x,4(i5,g12.4))
1767  WRITE(*,103) concut
1768 103 FORMAT(10x,' Cut on rms value is',g8.1)
1769  END IF
1770 
1771  IF(iact == 0) THEN
1772  WRITE(*,*) ' '
1773  WRITE(*,*) 'Improve constraints'
1774  END IF
1775  iact=1
1776 
1777  WRITE(*,102) iter,sum1,sum2,sum3
1778 102 FORMAT(i6,' rms',g12.4,' avrg_abs',g12.4,' max_abs',g12.4)
1779 
1780  CALL mpalloc(veccorrections,int(nvgb,mpl),'constraint corrections')
1781  veccorrections=0.0_mpd
1782 
1783  ! multiply (block-wise) inverse matrix and constraint vector
1784  isblck=0
1785  DO iblck=1,ncblck
1786  ifirst=matconsblocks(1,iblck) ! first constraint in block
1787  ilast=matconsblocks(1,iblck+1)-1 ! last constraint in block
1788  ncon=ilast+1-ifirst
1789  ieblck=isblck+(ncon*(ncon+1))/2
1790  CALL dbsvx(matconsproduct(isblck+1:ieblck),vecconsresiduals(ifirst:ilast),vecconssolution(ifirst:ilast),ncon)
1791  isblck=ieblck
1792  END DO
1793 
1794  DO jcgb=1,ncgb
1795  icgb=matconssort(3,jcgb) ! unsorted constraint index
1796  i=vecconsstart(icgb)
1797  rhs=listconstraints(i )%value ! right hand side
1798  sgm=listconstraints(i+1)%value ! sigma parameter
1799  DO j=i+2,vecconsstart(icgb+1)-1
1800  label=listconstraints(j)%label
1801  factr=listconstraints(j)%value
1802  itgbi=inone(label) ! -> ITGBI= index of parameter label
1803  ivgb =globalparlabelindex(2,itgbi) ! -> variable-parameter index
1804  IF(ivgb > 0) THEN
1805  veccorrections(ivgb)=veccorrections(ivgb)+vecconssolution(jcgb)*factr
1806  END IF
1807  ENDDO
1808  END DO
1809 
1810  DO i=1,nvgb ! add corrections
1811  itgbi=globalparvartototal(i)
1812  globalparameter(itgbi)=globalparameter(itgbi)+veccorrections(i)
1813  END DO
1814 
1815  CALL mpdealloc(veccorrections)
1816 
1817  END DO ! iteration 1 and 2
1818 
1819 END SUBROUTINE feasib ! make parameters feasible
1820 
1853 SUBROUTINE peread(more)
1854  USE mpmod
1855 
1856  IMPLICIT NONE
1857  INTEGER(mpi) :: i
1858  INTEGER(mpi) :: iact
1859  INTEGER(mpi) :: ierrc
1860  INTEGER(mpi) :: ierrf
1861  INTEGER(mpi) :: inder
1862  INTEGER(mpi) :: ioffp
1863  INTEGER(mpi) :: ios
1864  INTEGER(mpi) :: ithr
1865  INTEGER(mpi) :: jfile
1866  INTEGER(mpi) :: jrec
1867  INTEGER(mpi) :: k
1868  INTEGER(mpi) :: kfile
1869  INTEGER(mpi) :: l
1870  INTEGER(mpi) :: lun
1871  INTEGER(mpi) :: mpri
1872  INTEGER(mpi) :: n
1873  INTEGER(mpi) :: nact
1874  INTEGER(mpi) :: nbuf
1875  INTEGER(mpi) :: ndata
1876  INTEGER(mpi) :: noff
1877  INTEGER(mpi) :: noffs
1878  INTEGER(mpi) :: npointer
1879  INTEGER(mpi) :: npri
1880  INTEGER(mpi) :: nr
1881  INTEGER(mpi) :: nrc
1882  INTEGER(mpi) :: nrd
1883  INTEGER(mpi) :: nrpr
1884  INTEGER(mpi) :: nthr
1885  INTEGER(mpi) :: ntot
1886  INTEGER(mpi) :: maxRecordSize
1887  INTEGER(mpi) :: maxRecordFile
1888 
1889  INTEGER(mpi), INTENT(OUT) :: more
1890 
1891  LOGICAL :: lprint
1892  LOGICAL :: floop
1893  LOGICAL :: eof
1894  REAL(mpd) :: ds0
1895  REAL(mpd) :: ds1
1896  REAL(mpd) :: ds2
1897  REAL(mpd) :: dw
1898  !$ INTEGER(mpi) :: OMP_GET_THREAD_NUM
1899  CHARACTER (LEN=7) :: cfile
1900  SAVE
1901 
1902  inder(i)=readbufferdatai(i)
1903 
1904  DATA lprint/.true./
1905  DATA floop/.true./
1906  DATA npri / 0 /, mpri / 1000 /
1907  ! ...
1908  IF(ifile == 0) THEN ! start/restart
1909  nrec=0
1910  nrecd=0
1911  ntot=0
1912  sumrecords=0
1913  skippedrecords=0
1914  numblocks=0
1917  readbufferinfo=0 ! reset management info
1918  nrpr=1
1919  nthr=mthrdr
1920  nact=0 ! active threads (have something still to read)
1921  DO k=1,nthr
1922  IF (ifile < nfilb) THEN
1923  ifile=ifile+1
1924  readbufferinfo(1,k)=ifile
1925  readbufferinfo(2,k)=nact
1926  nact=nact+1
1927  END IF
1928  END DO
1929  END IF
1930  npointer=size(readbufferpointer)/nact
1931  ndata=size(readbufferdatai)/nact
1932  more=-1
1933  DO k=1,nthr
1934  iact=readbufferinfo(2,k)
1935  readbufferinfo(4,k)=0 ! reset counter
1936  readbufferinfo(5,k)=iact*ndata ! reset offset
1937  END DO
1938  numblocks=numblocks+1 ! new block
1939 
1940  !$OMP PARALLEL &
1941  !$OMP DEFAULT(PRIVATE) &
1942  !$OMP SHARED(readBufferInfo,readBufferPointer,readBufferDataI,readBufferDataD, &
1943  !$OMP readBufferDataF,nPointer,nData,skippedRecords,ndimbuf,NTHR,NFILF,FLOOP, &
1944  !$OMP IFD,KFD,IFILE,NFILB,WFD,XFD,icheck,keepOpen,ireeof,nrderr) &
1945  !$OMP NUM_THREADS(NTHR)
1946 
1947  ithr=1
1948  !$ ITHR=OMP_GET_THREAD_NUM()+1 ! thread number
1949  jfile=readbufferinfo(1,ithr) ! file index
1950  iact =readbufferinfo(2,ithr) ! active thread number
1951  jrec =readbufferinfo(3,ithr) ! records read
1952  ioffp=iact*npointer
1953  noffs=(ithr-1)*ndimbuf ! offset for intermediate float buffer
1954 
1955  files: DO WHILE (jfile > 0)
1956  kfile=kfd(2,jfile)
1957  ! open again
1958  IF (keepopen < 1 .AND. readbufferinfo(3,ithr) == 0) THEN
1959  CALL binopn(kfile,ithr,ios)
1960  END IF
1961  records: DO
1962  nbuf=readbufferinfo(4,ithr)+1
1963  noff=readbufferinfo(5,ithr)+2 ! 2 header words per record
1964  nr=ndimbuf
1965  IF(kfile <= nfilf) THEN ! Fortran file
1966  lun=kfile+10
1967  READ(lun,iostat=ierrf) n,(readbufferdataf(noffs+i),i=1,min(n/2,nr)),&
1968  (readbufferdatai(noff+i),i=1,min(n/2,nr))
1969  nr=n/2
1970  ! convert to double
1971  DO i=1,nr
1972  readbufferdatad(noff+i)=real(readbufferdataf(noffs+i),mpr8)
1973  END DO
1974  ! IF (ierrf < 0) REWIND lun ! end-of-file ! CHK use binrwd()
1975  eof=(ierrf /= 0)
1976  ELSE ! C file
1977  lun=kfile-nfilf
1978  IF (keepopen < 1) lun=ithr
1979 #ifdef READ_C_FILES
1980  CALL readc(readbufferdatad(noff+1),readbufferdataf(noffs+1),readbufferdatai(noff+1),nr,lun,ierrc)
1981  n=nr+nr
1982  IF (ierrc > 4) readbufferinfo(6,ithr)=readbufferinfo(6,ithr)+1
1983 #else
1984  ierrc=0
1985 #endif
1986  eof=(ierrc <= 0.AND.ierrc /= -4) ! allow buffer overruns -> skip record
1987  IF(eof.AND.ierrc < 0) THEN
1988  WRITE(*,*) 'Read error for binary Cfile', kfile, 'record', jrec+1, ':', ierrc
1989  WRITE(8,*) 'Read error for binary Cfile', kfile, 'record', jrec+1, ':', ierrc
1990  IF (icheck <= 0 .AND. ireeof <=0) THEN ! stop unless 'checkinput' mode or 'readerroraseof'
1991  WRITE(cfile,'(I7)') kfile
1992  CALL peend(18,'Aborted, read error(s) for binary file ' // cfile)
1993  stop 'PEREAD: stopping due to read errors'
1994  END IF
1995  IF (kfd(1,jfile) == 1) THEN ! count files with read errors in first loop
1996  !$OMP ATOMIC
1997  nrderr=nrderr+1
1998  END IF
1999  END IF
2000  END IF
2001  IF(eof) EXIT records ! end-of-files or error
2002 
2003  jrec=jrec+1
2004  readbufferinfo(3,ithr)=jrec
2005  IF(floop) THEN
2006  xfd(jfile)=max(xfd(jfile),n)
2007  IF(ithr == 1) THEN
2008  CALL hmplnt(1,n)
2009  IF(inder(noff+1) /= 0) CALL hmpent(8,real(inder(noff+1),mps))
2010  END IF
2011  END IF
2012 
2013  IF (nr <= ndimbuf) THEN
2014  readbufferinfo(4,ithr)=nbuf
2015  readbufferinfo(5,ithr)=noff+nr
2016 
2017  readbufferpointer(ioffp+nbuf)=noff ! pointer to start of buffer
2018  readbufferdatai(noff )=noff+nr ! pointer to end of buffer
2019  readbufferdatai(noff-1)=ifd(kfile)+jrec ! global record number (available with LOOP2)
2020  readbufferdatad(noff )=real(kfile,mpr8) ! file number
2021  readbufferdatad(noff-1)=real(wfd(kfile),mpr8) ! weight
2022 
2023  IF ((noff+nr+2+ndimbuf >= ndata*(iact+1)).OR.(nbuf >= npointer)) EXIT files ! buffer full
2024  ELSE
2025  !$OMP ATOMIC
2027  cycle records
2028  END IF
2029 
2030  END DO records
2031 
2032  readbufferinfo(1,ithr)=-jfile ! flag eof
2033  IF (keepopen < 1) THEN ! close again
2034  CALL bincls(kfile,ithr)
2035  ELSE ! rewind
2036  CALL binrwd(kfile)
2037  END IF
2038  IF (kfd(1,jfile) == 1) THEN
2039  print *, 'PEREAD: file ', kfile, 'read the first time, found',jrec,' records'
2040  kfd(1,jfile)=-jrec
2041  ELSE
2042  !PRINT *, 'PEREAD: file ', kfile, 'records', jrec, -kfd(1,jfile)
2043  IF (-kfd(1,jfile) /= jrec) THEN
2044  WRITE(cfile,'(I7)') kfile
2045  CALL peend(19,'Aborted, binary file modified (length) ' // cfile)
2046  stop 'PEREAD: file modified (length)'
2047  END IF
2048  END IF
2049  ! take next file
2050  !$OMP CRITICAL
2051  IF (ifile < nfilb) THEN
2052  ifile=ifile+1
2053  jrec=0
2054  readbufferinfo(1,ithr)=ifile
2055  readbufferinfo(3,ithr)=jrec
2056  END IF
2057  !$OMP END CRITICAL
2058  jfile=readbufferinfo(1,ithr)
2059 
2060  END DO files
2061  !$OMP END PARALLEL
2062  ! compress pointers
2063  nrd=readbufferinfo(4,1) ! buffers from 1 .thread
2064  DO k=2,nthr
2065  iact =readbufferinfo(2,k)
2066  ioffp=iact*npointer
2067  nbuf=readbufferinfo(4,k)
2068  DO l=1,nbuf
2069  readbufferpointer(nrd+l)=readbufferpointer(ioffp+l)
2070  END DO
2071  nrd=nrd+nbuf
2072  END DO
2073 
2074  more=0
2075  DO k=1,nthr
2076  jfile=readbufferinfo(1,k)
2077  IF (jfile > 0) THEN ! no eof yet
2078  readbufferinfo(2,k)=more
2079  more=more+1
2080  ELSE
2081  ! no more files, thread retires
2082  readbufferinfo(1,k)=0
2083  readbufferinfo(2,k)=-1
2084  readbufferinfo(3,k)=0
2086  readbufferinfo(6,k)=0
2087  END IF
2088  END DO
2089  ! record limit ?
2090  IF (mxrec > 0.AND.(ntot+nrd) >= mxrec) THEN
2091  nrd=mxrec-ntot
2092  more=-1
2093  DO k=1,nthr
2094  jfile=readbufferinfo(1,k)
2095  IF (jfile > 0) THEN ! rewind or close files
2096  nrc=readbufferinfo(3,k)
2097  IF (kfd(1,jfile) == 1) kfd(1,jfile)=-nrc
2098  kfile=kfd(2,jfile)
2099  IF (keepopen < 1) THEN ! close again
2100  CALL bincls(kfile,k)
2101  ELSE ! rewind
2102  CALL binrwd(kfile)
2103  END IF
2104  END IF
2105  END DO
2106  END IF
2107 
2108  ntot=ntot+nrd
2109  nrec=ntot
2110  numreadbuffer=nrd
2111 
2115 
2116  DO WHILE (nloopn == 0.AND.ntot >= nrpr)
2117  WRITE(*,*) ' Record ',nrpr
2118  IF (nrpr < 100000) THEN
2119  nrpr=nrpr*10
2120  ELSE
2121  nrpr=nrpr+100000
2122  END IF
2123  END DO
2124 
2125  IF (ncache > 0.AND.nloopn <= 1.AND. npri < mpri.AND.mprint > 1) THEN
2126  npri=npri+1
2127  IF (npri == 1) WRITE(*,100)
2128  WRITE(*,101) nrec, nrd, more ,ifile
2129 100 FORMAT(/' PeRead records active file' &
2130  /' total block threads number')
2131 101 FORMAT(' PeRead',4i10)
2132  END IF
2133 
2134  IF (more <= 0) THEN
2135  ifile=0
2136  IF (floop) THEN
2137  ! check for file weights
2138  ds0=0.0_mpd
2139  ds1=0.0_mpd
2140  ds2=0.0_mpd
2141  maxrecordsize=0
2142  maxrecordfile=0
2143  DO k=1,nfilb
2144  IF (xfd(k) > maxrecordsize) THEN
2145  maxrecordsize=xfd(k)
2146  maxrecordfile=k
2147  END IF
2148  dw=real(-kfd(1,k),mpd)
2149  IF (wfd(k) /= 1.0) nfilw=nfilw+1
2150  ds0=ds0+dw
2151  ds1=ds1+dw*real(wfd(k),mpd)
2152  ds2=ds2+dw*real(wfd(k)**2,mpd)
2153  END DO
2154  print *, 'PEREAD: file ', maxrecordfile, 'with max record size ', maxrecordsize
2155  IF (nfilw > 0.AND.ds0 > 0.0_mpd) THEN
2156  ds1=ds1/ds0
2157  ds2=ds2/ds0-ds1*ds1
2158  DO lun=6,lunlog,2
2159  WRITE(lun,177) nfilw,real(ds1,mps),real(ds2,mps)
2160 177 FORMAT(/' !!!!!',i4,' weighted binary files', &
2161  /' !!!!! mean, variance of weights =',2g12.4)
2162  END DO
2163  END IF
2164  ! integrate record numbers
2165  DO k=2,nfilb
2166  ifd(k)=ifd(k-1)-kfd(1,k-1)
2167  END DO
2168  ! sort
2169  IF (nthr > 1) CALL sort2k(kfd,nfilb)
2170  IF (skippedrecords > 0) THEN
2171  print *, 'PEREAD skipped records: ', skippedrecords
2172  ndimbuf=maxrecordsize/2 ! adjust buffer size
2173  END IF
2174  END IF
2175  lprint=.false.
2176  floop=.false.
2177  IF (ncache > 0.AND.nloopn <= 1.AND.mprint > 0) &
2179 179 FORMAT(/' Read cache usage (#blocks, #records, ', &
2180  'min,max records/block'/17x,4i10)
2181  END IF
2182  RETURN
2183 
2184 END SUBROUTINE peread
2185 
2193 SUBROUTINE peprep(mode)
2194  USE mpmod
2195 
2196  IMPLICIT NONE
2197 
2198  INTEGER(mpi), INTENT(IN) :: mode
2199 
2200  INTEGER(mpi) :: ibuf
2201  INTEGER(mpi) :: ichunk
2202  INTEGER(mpi) :: isfrst
2203  INTEGER(mpi) :: islast
2204  INTEGER(mpi) :: ist
2205  INTEGER(mpi) :: itgbi
2206  INTEGER(mpi) :: j
2207  INTEGER(mpi) :: ja
2208  INTEGER(mpi) :: jb
2209  INTEGER(mpi) :: jsp
2210  INTEGER(mpi) :: nst
2211  INTEGER(mpi), PARAMETER :: maxbad = 100 ! max number of bad records with print out
2212  INTEGER(mpi) :: nbad
2213  INTEGER(mpi) :: nerr
2214  INTEGER(mpi) :: inone
2215 
2216 
2217  isfrst(ibuf)=readbufferpointer(ibuf)+1
2218  islast(ibuf)=readbufferdatai(readbufferpointer(ibuf))
2219 
2220  IF (mode > 0) THEN
2221 #ifdef __PGIC__
2222  ! to prevent "PGF90-F-0000-Internal compiler error. Could not locate uplevel instance for stblock"
2223  ichunk=256
2224 #else
2225  ichunk=min((numreadbuffer+mthrd-1)/mthrd/32+1,256)
2226 #endif
2227  ! parallelize record loop
2228  !$OMP PARALLEL DO &
2229  !$OMP DEFAULT(PRIVATE) &
2230  !$OMP SHARED(numReadBuffer,readBufferPointer,readBufferDataI,readBufferDataD,ICHUNK,iscerr,dscerr) &
2231  !$OMP SCHEDULE(DYNAMIC,ICHUNK)
2232  DO ibuf=1,numreadbuffer ! buffer for current record
2233  ist=isfrst(ibuf)
2234  nst=islast(ibuf)
2235  DO ! loop over measurements
2236  CALL isjajb(nst,ist,ja,jb,jsp)
2237  IF(jb == 0) EXIT
2238  DO j=1,ist-jb
2239  readbufferdatai(jb+j)=inone( readbufferdatai(jb+j) ) ! translate to index
2240  END DO
2241  ! scale error ?
2242  IF (iscerr > 0) THEN
2243  IF (jb < ist) THEN
2244  readbufferdatad(jb) = readbufferdatad(jb) * dscerr(1) ! 'global' measurement
2245  ELSE
2246  readbufferdatad(jb) = readbufferdatad(jb) * dscerr(2) ! 'local' measurement
2247  END IF
2248  END IF
2249  END DO
2250  END DO
2251  !$OMP END PARALLEL DO
2252  END IF
2253 
2254  !$POMP INST BEGIN(peprep)
2255  IF (mode <= 0) THEN
2256  nbad=0
2257  DO ibuf=1,numreadbuffer ! buffer for current record
2258  CALL pechk(ibuf,nerr)
2259  IF(nerr > 0) THEN
2260  nbad=nbad+1
2261  IF(nbad >= maxbad) EXIT
2262  ELSE
2263  ist=isfrst(ibuf)
2264  nst=islast(ibuf)
2265  DO ! loop over measurements
2266  CALL isjajb(nst,ist,ja,jb,jsp)
2267  IF(jb == 0) EXIT
2268  DO j=1,ist-jb
2269  itgbi=inone( readbufferdatai(jb+j) ) ! generate index
2270  END DO
2271  END DO
2272  END IF
2273  END DO
2274  IF(nbad > 0) THEN
2275  CALL peend(20,'Aborted, bad binary records')
2276  stop 'PEREAD: stopping due to bad records'
2277  END IF
2278  END IF
2279  !$POMP INST END(peprep)
2280 
2281 END SUBROUTINE peprep
2282 
2290 SUBROUTINE pechk(ibuf, nerr)
2291  USE mpmod
2292 
2293  IMPLICIT NONE
2294  REAL(mpr8) :: glder
2295  INTEGER(mpi) :: i
2296  INTEGER(mpi) :: is
2297  INTEGER(mpi) :: ist
2298  INTEGER(mpi) :: inder
2299  INTEGER(mpi) :: ioff
2300  INTEGER(mpi) :: isfrst
2301  INTEGER(mpi) :: islast
2302  INTEGER(mpi) :: ja
2303  INTEGER(mpi) :: jb
2304  INTEGER(mpi) :: jsp
2305  INTEGER(mpi) :: nan
2306  INTEGER(mpi) :: nst
2307 
2308  INTEGER(mpi), INTENT(IN) :: ibuf
2309  INTEGER(mpi), INTENT(OUT) :: nerr
2310  SAVE
2311  ! ...
2312  inder(i)=readbufferdatai(i)
2313  glder(i)=readbufferdatad(i)
2314  isfrst(ibuf)=readbufferpointer(ibuf)+1
2315  islast(ibuf)=readbufferdatai(readbufferpointer(ibuf))
2316 
2317  ist=isfrst(ibuf)
2318  nst=islast(ibuf)
2319  nerr=0
2320  is=ist
2321  jsp=0
2322  outer: DO WHILE(is < nst)
2323  ja=0
2324  jb=0
2325  inner1: DO
2326  is=is+1
2327  IF(is > nst) EXIT outer
2328  IF(inder(is) == 0) EXIT inner1 ! found 1. marker
2329  END DO inner1
2330  ja=is
2331  inner2: DO
2332  is=is+1
2333  IF(is > nst) EXIT outer
2334  IF(inder(is) == 0) EXIT inner2 ! found 2. marker
2335  END DO inner2
2336  jb=is
2337  IF(ja+1 == jb.AND.glder(jb) < 0.0_mpr8) THEN
2338  ! special data
2339  jsp=jb ! pointer to special data
2340  is=is+nint(-glder(jb),mpi) ! skip NSP words
2341  cycle outer
2342  END IF
2343  DO WHILE(inder(is+1) /= 0.AND.is < nst)
2344  is=is+1
2345  END DO
2346  END DO outer
2347  IF(is > nst) THEN
2348  ioff = readbufferpointer(ibuf)
2349  WRITE(*,100) readbufferdatai(ioff-1), int(readbufferdatad(ioff),mpi)
2350 100 FORMAT(' PEREAD: record ', i8,' in file ',i6, ' is broken !!!')
2351  nerr=nerr+1
2352  ENDIF
2353  nan=0
2354  DO i=ist, nst
2355  IF(.NOT.(readbufferdatad(i) <= 0.0_mpr8).AND..NOT.(readbufferdatad(i) > 0.0_mpr8)) nan=nan+1
2356  END DO
2357  IF(nan > 0) THEN
2358  ioff = readbufferpointer(ibuf)
2359  WRITE(*,101) readbufferdatai(ioff-1), int(readbufferdatad(ioff),mpi), nan
2360 101 FORMAT(' PEREAD: record ', i8,' in file ',i6, ' contains ', i6, ' NaNs !!!')
2361  nerr= nerr+2
2362  ENDIF
2363 
2364 END SUBROUTINE pechk
2365 
2370 SUBROUTINE pepgrp
2371  USE mpmod
2372  USE mpdalc
2373 
2374  IMPLICIT NONE
2375 
2376  INTEGER(mpi) :: ibuf
2377  INTEGER(mpi) :: ichunk
2378  INTEGER(mpi) :: iproc
2379  INTEGER(mpi) :: ioff
2380  INTEGER(mpi) :: ioffbi
2381  INTEGER(mpi) :: isfrst
2382  INTEGER(mpi) :: islast
2383  INTEGER(mpi) :: ist
2384  INTEGER(mpi) :: itgbi
2385  INTEGER(mpi) :: j
2386  INTEGER(mpi) :: ja
2387  INTEGER(mpi) :: jb
2388  INTEGER(mpi) :: jsp
2389  INTEGER(mpi) :: nalg
2390  INTEGER(mpi) :: nst
2391  INTEGER(mpi) :: inone
2392  INTEGER(mpl) :: length
2393  !$ INTEGER(mpi) :: OMP_GET_THREAD_NUM
2394 
2395  isfrst(ibuf)=readbufferpointer(ibuf)+1
2396  islast(ibuf)=readbufferdatai(readbufferpointer(ibuf))
2397 
2398  CALL useone ! make (INONE) usable
2399  globalparheader(-2)=-1 ! set flag to inhibit further updates
2400  ! need back index
2401  IF (mcount > 0) THEN
2402  length=globalparheader(-1)*mthrd
2403  CALL mpalloc(backindexusage,length,'global variable-index array')
2404  backindexusage=0
2405  END IF
2406 
2407 #ifdef __PGIC__
2408  ! to prevent "PGF90-F-0000-Internal compiler error. Could not locate uplevel instance for stblock"
2409  ichunk=256
2410 #else
2411  ichunk=min((numreadbuffer+mthrd-1)/mthrd/32+1,256)
2412 #endif
2413  ! parallelize record loop
2414  !$OMP PARALLEL DO &
2415  !$OMP DEFAULT(PRIVATE) &
2416  !$OMP SHARED(numReadBuffer,readBufferPointer,readBufferDataI,readBufferDataD,backIndexUsage,globalParHeader,ICHUNK,MCOUNT) &
2417  !$OMP SCHEDULE(DYNAMIC,ICHUNK)
2418  DO ibuf=1,numreadbuffer ! buffer for current record
2419  ist=isfrst(ibuf)
2420  nst=islast(ibuf)
2421  IF (mcount > 0) THEN
2422  ! count per record
2423  iproc=0
2424  !$ IPROC=OMP_GET_THREAD_NUM() ! thread number
2425  ioffbi=globalparheader(-1)*iproc
2426  nalg=0
2427  ioff=readbufferpointer(ibuf)
2428  DO ! loop over measurements
2429  CALL isjajb(nst,ist,ja,jb,jsp)
2430  IF(jb == 0) EXIT
2431  IF (ist > jb) THEN
2432  DO j=1,ist-jb
2433  itgbi=inone( readbufferdatai(jb+j) ) ! translate to index
2434  IF (backindexusage(ioffbi+itgbi) == 0) THEN
2435  nalg=nalg+1
2436  readbufferdatai(ioff+nalg)=itgbi
2437  backindexusage(ioffbi+itgbi)=nalg
2438  END IF
2439  END DO
2440  END IF
2441  END DO
2442  ! reset back index
2443  DO j=1,nalg
2444  itgbi=readbufferdatai(ioff+j)
2445  backindexusage(ioffbi+itgbi)=0
2446  END DO
2447  ! sort (record)
2448  CALL sort1k(readbufferdatai(ioff+1),nalg)
2449  readbufferdatai(ioff)=ioff+nalg
2450  ELSE
2451  ! count per equation
2452  DO ! loop over measurements
2453  CALL isjajb(nst,ist,ja,jb,jsp)
2454  IF(jb == 0) EXIT
2455  IF (ist > jb) THEN
2456  DO j=1,ist-jb
2457  readbufferdatai(jb+j)=inone( readbufferdatai(jb+j) ) ! translate to index
2458  END DO
2459  ! sort (equation)
2460  CALL sort1k(readbufferdatai(jb+1),ist-jb)
2461  END IF
2462  END DO
2463  END IF
2464  END DO
2465  !$OMP END PARALLEL DO
2466 
2467  !$POMP INST BEGIN(pepgrp)
2468  DO ibuf=1,numreadbuffer ! buffer for current record
2469  ist=isfrst(ibuf)
2470  nst=islast(ibuf)
2471  IF (mcount == 0) THEN
2472  ! equation level
2473  DO ! loop over measurements
2474  CALL isjajb(nst,ist,ja,jb,jsp)
2475  IF(jb == 0) EXIT
2476  CALL pargrp(jb+1,ist)
2477  END DO
2478  ELSE
2479  ! record level, group
2480  CALL pargrp(ist,nst)
2481  ! count
2482  DO j=ist,nst
2483  itgbi=readbufferdatai(j)
2484  globalparlabelindex(4,itgbi)=globalparlabelindex(4,itgbi)+1
2485  END DO
2486  ENDIF
2487  END DO
2488  ! free back index
2489  IF (mcount > 0) THEN
2491  END IF
2492  !$POMP INST END(pepgrp)
2493  globalparheader(-2)=0 ! reset flag to reenable further updates
2494 
2495 END SUBROUTINE pepgrp
2496 
2504 SUBROUTINE pargrp(inds,inde)
2505  USE mpmod
2506 
2507  IMPLICIT NONE
2508 
2509  INTEGER(mpi) :: istart
2510  INTEGER(mpi) :: itgbi
2511  INTEGER(mpi) :: j
2512  INTEGER(mpi) :: jstart
2513  INTEGER(mpi) :: jtgbi
2514  INTEGER(mpi) :: lstart
2515  INTEGER(mpi) :: ltgbi
2516 
2517  INTEGER(mpi), INTENT(IN) :: inds
2518  INTEGER(mpi), INTENT(IN) :: inde
2519 
2520  IF (inds > inde) RETURN
2521  ltgbi=-1
2522  lstart=-1
2523  ! build up groups
2524  DO j=inds,inde
2525  itgbi=readbufferdatai(j)
2526  istart=globalparlabelindex(3,itgbi) ! label of group start
2527  IF (istart == 0) THEN ! not yet in group
2528  IF (itgbi /= ltgbi+1)THEN ! start group
2529  globalparlabelindex(3,itgbi)=globalparlabelindex(1,itgbi)
2530  ELSE ! extend group
2531  globalparlabelindex(3,itgbi)=globalparlabelindex(3,ltgbi)
2532  END IF
2533  END IF
2534  ltgbi=itgbi
2535  END DO
2536  ! split groups:
2537  ! - start inside group?
2538  itgbi=readbufferdatai(inds)
2539  istart=globalparlabelindex(3,itgbi) ! label of group start
2540  jstart=globalparlabelindex(1,itgbi) ! label of first parameter
2541  IF (istart /= jstart) THEN ! start new group
2542  DO WHILE (globalparlabelindex(3,itgbi) == istart)
2543  globalparlabelindex(3,itgbi) = jstart
2544  itgbi=itgbi+1
2545  IF (itgbi > globalparheader(-1)) EXIT
2546  END DO
2547  END IF
2548  ! - not neigbours anymore
2549  ltgbi=readbufferdatai(inds)
2550  DO j=inds+1,inde
2551  itgbi=readbufferdatai(j)
2552  IF (itgbi /= ltgbi+1) THEN
2553  ! split after ltgbi
2554  lstart=globalparlabelindex(3,ltgbi) ! label of last group start
2555  jtgbi=ltgbi+1 ! new group after ltgbi
2556  jstart=globalparlabelindex(1,jtgbi)
2557  DO WHILE (globalparlabelindex(3,jtgbi) == lstart)
2558  globalparlabelindex(3,jtgbi) = jstart
2559  jtgbi=jtgbi+1
2560  IF (jtgbi > globalparheader(-1)) EXIT
2561  IF (jtgbi == itgbi) jstart=globalparlabelindex(1,jtgbi)
2562  END DO
2563  ! split at itgbi
2564  jtgbi=itgbi
2565  istart=globalparlabelindex(3,jtgbi) ! label of group start
2566  jstart=globalparlabelindex(1,jtgbi) ! label of first parameter
2567  IF (istart /= jstart) THEN ! start new group
2568  DO WHILE (globalparlabelindex(3,jtgbi) == istart)
2569  globalparlabelindex(3,jtgbi) = jstart
2570  jtgbi=jtgbi+1
2571  IF (jtgbi > globalparheader(-1)) EXIT
2572  END DO
2573  END IF
2574  ENDIF
2575  ltgbi=itgbi
2576  END DO
2577  ! - end inside group?
2578  itgbi=readbufferdatai(inde)
2579  IF (itgbi < globalparheader(-1)) THEN
2580  istart=globalparlabelindex(3,itgbi) ! label of group start
2581  itgbi=itgbi+1
2582  jstart=globalparlabelindex(1,itgbi) ! label of new group start
2583  DO WHILE (globalparlabelindex(3,itgbi) == istart)
2584  globalparlabelindex(3,itgbi) = jstart
2585  itgbi=itgbi+1
2586  IF (itgbi > globalparheader(-1)) EXIT
2587  END DO
2588  END IF
2589 
2590 END SUBROUTINE pargrp
2591 
2614 SUBROUTINE isjajb(nst,is,ja,jb,jsp)
2615  USE mpmod
2616 
2617  IMPLICIT NONE
2618  REAL(mpr8) :: glder
2619  INTEGER(mpi) :: i
2620  INTEGER(mpi) :: inder
2621 
2622  INTEGER(mpi), INTENT(IN) :: nst
2623  INTEGER(mpi), INTENT(IN OUT) :: is
2624  INTEGER(mpi), INTENT(OUT) :: ja
2625  INTEGER(mpi), INTENT(OUT) :: jb
2626  INTEGER(mpi), INTENT(OUT) :: jsp
2627  SAVE
2628  ! ...
2629  inder(i)=readbufferdatai(i)
2630  glder(i)=readbufferdatad(i)
2631 
2632  jsp=0
2633  DO
2634  ja=0
2635  jb=0
2636  IF(is >= nst) RETURN
2637  DO
2638  is=is+1
2639  IF(inder(is) == 0) EXIT
2640  END DO
2641  ja=is
2642  DO
2643  is=is+1
2644  IF(inder(is) == 0) EXIT
2645  END DO
2646  jb=is
2647  IF(ja+1 == jb.AND.glder(jb) < 0.0_mpr8) THEN
2648  ! special data
2649  jsp=jb ! pointer to special data
2650  is=is+nint(-glder(jb),mpi) ! skip NSP words
2651  cycle
2652  END IF
2653  DO WHILE(inder(is+1) /= 0.AND.is < nst)
2654  is=is+1
2655  END DO
2656  EXIT
2657  END DO
2658 
2659 END SUBROUTINE isjajb
2660 
2661 
2662 !***********************************************************************
2663 ! LOOPN ...
2669 
2670 SUBROUTINE loopn
2671  USE mpmod
2672 
2673  IMPLICIT NONE
2674  REAL(mpd) :: dsum
2675  REAL(mps) :: elmt
2676  REAL(mpd) :: factrj
2677  REAL(mpd) :: factrk
2678  REAL(mpr8) :: glder
2679  REAL(mps) :: peakd
2680  REAL(mps) :: peaki
2681  REAL(mps) :: ratae
2682  REAL(mpd) :: rhs
2683  REAL(mps) :: rloop
2684  REAL(mpd) :: sgm
2685  REAL(mps) :: used
2686  REAL(mps) :: usei
2687  REAL(mpd) :: weight
2688  INTEGER(mpi) :: i
2689  INTEGER(mpi) :: ia
2690  INTEGER(mpi) :: ib
2691  INTEGER(mpi) :: ibuf
2692  INTEGER(mpi) :: inder
2693  INTEGER(mpi) :: ioffb
2694  INTEGER(mpi) :: ipr
2695  INTEGER(mpi) :: isfrst
2696  INTEGER(mpi) :: islast
2697  INTEGER(mpi) :: itgbi
2698  INTEGER(mpi) :: itgbij
2699  INTEGER(mpi) :: itgbik
2700  INTEGER(mpi) :: ivgb
2701  INTEGER(mpi) :: ivgbij
2702  INTEGER(mpi) :: ivgbik
2703  INTEGER(mpi) :: j
2704  INTEGER(mpi) :: k
2705  INTEGER(mpi) :: lastit
2706  INTEGER(mpi) :: lun
2707  INTEGER(mpi) :: ncrit
2708  INTEGER(mpi) :: ndfs
2709  INTEGER(mpi) :: ngras
2710  INTEGER(mpi) :: nparl
2711  INTEGER(mpi) :: nr
2712  INTEGER(mpi) :: nrej
2713  INTEGER(mpi) :: inone
2714  INTEGER(mpi) :: ilow
2715  INTEGER(mpi) :: nlow
2716  INTEGER(mpi) :: nzero
2717  LOGICAL :: btest
2718 
2719  REAL(mpd):: adder
2720  REAL(mpd)::funref
2721  REAL(mpd)::dchi2s
2722  REAL(mpd)::matij
2723  REAL(mpd)::sndf
2724  SAVE
2725  ! ...
2726  isfrst(ibuf)=readbufferpointer(ibuf)+1
2727  islast(ibuf)=readbufferdatai(readbufferpointer(ibuf))
2728  inder(i)=readbufferdatai(i)
2729  glder(i)=readbufferdatad(i)
2730  ! ----- book and reset ---------------------------------------------
2731  IF(nloopn == 0) THEN ! first call
2732  lastit=-1
2733  iitera=0
2734  END IF
2735 
2736  nloopn=nloopn+1 ! increase loop counter
2737  ndfsum=0
2738  sumndf=0.0_mpd
2739  funref=0.0_mpd
2740 
2741  IF(nloopn == 1) THEN ! book histograms for 1. iteration
2742  CALL gmpdef(1,4,'Function value in iterations')
2743  IF (metsol == 4 .OR. metsol == 5) THEN ! extend to GMRES, i.e. 6?
2744  CALL gmpdef(2,3,'Number of MINRES steps vs iteration nr')
2745  END IF
2746  CALL hmpdef( 5,0.0,0.0,'Number of degrees of freedom')
2747  CALL hmpdef(11,0.0,0.0,'Number of local parameters')
2748  CALL hmpdef(23,0.0,0.0, 'SQRT of diagonal elements without presigma')
2749  CALL hmpdef(24,0.0,0.0, 'Log10 of off-diagonal elements')
2750  CALL hmpdef(25,0.0,0.0, 'Relative individual pre-sigma')
2751  CALL hmpdef(26,0.0,0.0, 'Relative global pre-sigma')
2752  END IF
2753 
2754 
2755  CALL hmpdef(3,-prange,prange, & ! book
2756  'Normalized residuals of single (global) measurement')
2757  CALL hmpdef(12,-prange,prange, & ! book
2758  'Normalized residuals of single (local) measurement')
2759  CALL hmpdef(13,-prange,prange, & ! book
2760  'Pulls of single (global) measurement')
2761  CALL hmpdef(14,-prange,prange, & ! book
2762  'Pulls of single (local) measurement')
2763  CALL hmpdef(4,0.0,0.0,'Chi^2/Ndf after local fit')
2764  CALL gmpdef(4,5,'location, dispersion (res.) vs record nr')
2765  CALL gmpdef(5,5,'location, dispersion (pull) vs record nr')
2766 
2767  ! WRITE(*,*) 'LOOPN ', NLOOPN, ' executing ICALCM=', ICALCM
2768 
2769  ! reset
2770 
2771  globalvector=0.0_mpd ! reset rhs vector IGVEC
2772  globalcounter=0
2773  IF(icalcm == 1) THEN
2774  globalmatd=0.0_mpd
2775  globalmatf=0.
2776  IF (metsol >= 4.AND.metsol < 7) matprecond=0.0_mpd
2777  END IF
2778 
2779  IF(nloopn == 2) CALL hmpdef(6,0.0,0.0,'Down-weight fraction')
2780 
2781  newite=.false.
2782  IF(iterat /= lastit) THEN ! new iteration
2783  newite=.true.
2784  funref=fvalue
2785  IF(nloopn > 1) THEN
2786  nrej=nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3)
2787  ! CALL MEND
2788  IF(iterat == 1) THEN
2789  chicut=chirem
2790  ELSE IF(iterat >= 1) THEN
2791  chicut=sqrt(chicut)
2792  IF(chicut /= 0.0.AND.chicut < 1.5) chicut=1.0
2793  IF(chicut /= 0.0.AND.nrej == 0) chicut=1.0
2794  END IF
2795  END IF
2796  ! WRITE(*,111) ! header line
2797  END IF
2798 
2799  DO i=0,3
2800  nrejec(i)=0 ! reset reject counter
2801  END DO
2802  DO k=3,6
2803  writebufferheader(k)=0 ! cache usage
2804  writebufferheader(-k)=0
2805  END DO
2806  ! statistics per binary file
2807  DO i=1,nfilb
2808  jfd(i)=0
2809  cfd(i)=0.0
2810  dfd(i)=0
2811  END DO
2812 
2813  IF (imonit /= 0) meashists=0 ! reset monitoring histograms
2814 
2815  ! ----- read next data ----------------------------------------------
2816  DO
2817  CALL peread(nr) ! read records
2818  CALL peprep(1) ! prepare records
2819  ndfs =0
2820  sndf =0.0_mpd
2821  dchi2s=0.0_mpd
2822  CALL loopbf(nrejec,ndfs,sndf,dchi2s,nfiles,jfd,cfd,dfd)
2823  ndfsum=ndfsum+ndfs
2824  sumndf=sumndf+sndf
2825  CALL addsum(dchi2s)
2826  IF (nr <= 0) EXIT ! next block of events ?
2827  END DO
2828  ! sum up RHS (over threads) once (reduction in LOOPBF: summation for each block)
2829  ioffb=0
2830  DO ipr=2,mthrd
2831  ioffb=ioffb+lenglobalvec
2832  DO k=1,lenglobalvec
2833  globalvector(k)=globalvector(k)+globalvector(ioffb+k)
2835  END DO
2836  END DO
2837 
2838  IF (icalcm == 1) THEN
2839  ! PRINT *, ' cache/w ',(writeBufferHeader(-K),K=3,6),(writeBufferHeader(K),K=3,6)
2840  nparl=writebufferheader(3)
2841  ncrit=writebufferheader(4)
2842  used=real(writebufferheader(-5),mps)/real(writebufferheader(-3),mps)*0.1
2843  usei=real(writebufferheader(5),mps)/real(writebufferheader(3),mps)*0.1
2844  peakd=real(writebufferheader(-6),mps)*0.1
2845  peaki=real(writebufferheader(6),mps)*0.1
2846  WRITE(*,111) nparl,ncrit,usei,used,peaki,peakd
2847 111 FORMAT(' Write cache usage (#flush,#overrun,<levels>,', &
2848  'peak(levels))'/2i7,',',4(f6.1,'%'))
2849  ! fill part of MINRES preconditioner matrix from binary files (formerly in mgupdt)
2850  IF (metsol >= 4.AND.metsol < 7) THEN
2851  IF (mbandw == 0) THEN
2852  ! default preconditioner (diagonal)
2853  DO i=1, nvgb
2854  matprecond(i)=matij(i,i)
2855  END DO
2856  ELSE IF (mbandw > 0) THEN
2857  ! band matrix
2858  DO i=1, nvgb
2859  ia=indprecond(i) ! index of diagonal element
2860  DO j=max(1,i-mbandw+1),i
2861  matprecond(ia-i+j)=matij(i,j)
2862  END DO
2863  END DO
2864  END IF
2865  END IF
2866  ! fill second half (j>i) of global matric for extended storage
2867  IF (mextnd > 0) CALL mhalf2()
2868  END IF
2869 
2870  ! check entries/counters
2871  nlow=0
2872  ilow=1
2873  nzero=0
2874  DO i=1,nvgb
2875  IF(globalcounter(i) == 0) nzero=nzero+1
2876  IF(globalcounter(i) < mreqena) THEN
2877  nlow=nlow+1
2878  IF(globalcounter(i) < globalcounter(ilow)) ilow=i
2879  END IF
2880  END DO
2881  IF(nlow > 0) THEN
2882  nalow=nalow+nlow
2883  IF(icalcm == 1) nxlow=max(nxlow,nlow) ! for matrix construction ?
2884  itgbi=globalparvartototal(ilow)
2885  print *
2886  print *, " ... warning ..."
2887  print *, " global parameters with too few (< MREQENA) accepted entries: ", nlow
2888  print *, " minimum entries: ", globalcounter(ilow), " for label ", globalparlabelindex(1,itgbi)
2889  print *
2890  END IF
2891  IF(icalcm == 1 .AND. nzero > 0) THEN
2892  ndefec = nzero ! rank defect
2893  WRITE(*,*) 'Warning: the rank defect of the symmetric',nfgb, &
2894  '-by-',nfgb,' matrix is ',ndefec,' (should be zero).'
2895  WRITE(lun,*) 'Warning: the rank defect of the symmetric',nfgb, &
2896  '-by-',nfgb,' matrix is ',ndefec,' (should be zero).'
2897  IF (iforce == 0) THEN
2898  isubit=1
2899  WRITE(*,*) ' --> enforcing SUBITO mode'
2900  WRITE(lun,*) ' --> enforcing SUBITO mode'
2901  END IF
2902  END IF
2903 
2904  ! ----- after end-of-data add contributions from pre-sigma ---------
2905 
2906  IF(nloopn == 1) THEN
2907  ! plot diagonal elements
2908  elmt=0.0
2909  DO i=1,nvgb ! diagonal elements
2910  elmt=real(matij(i,i),mps)
2911  IF(elmt > 0.0) CALL hmpent(23,1.0/sqrt(elmt))
2912  END DO
2913  END IF
2914 
2915 
2916 
2917  ! add pre-sigma contributions to matrix diagonal
2918 
2919  ! WRITE(*,*) 'Adding to diagonal ICALCM IND6',ICALCM,IND6
2920 
2921  IF(icalcm == 1) THEN
2922  DO ivgb=1,nvgb ! add evtl. pre-sigma
2923  ! WRITE(*,*) 'Index ',IVGB,IVGB,QM(IND6+IVGB)
2924  IF(globalparpreweight(ivgb) /= 0.0) THEN
2925  IF(ivgb > 0) CALL mupdat(ivgb,ivgb,globalparpreweight(ivgb))
2926  END IF
2927  END DO
2928  END IF
2929 
2930  CALL hmpwrt(23)
2931  CALL hmpwrt(24)
2932  CALL hmpwrt(25)
2933  CALL hmpwrt(26)
2934 
2935 
2936  ! add regularization term to F and to rhs --------------------------
2937 
2938  ! WRITE(*,*) 'NREGUL ',NREGUL,NLOOPN
2939 
2940  IF(nregul /= 0) THEN ! add regularization term to F and to rhs
2941  DO ivgb=1,nvgb
2942  itgbi=globalparvartototal(ivgb) ! global parameter index
2944  adder=globalparpreweight(ivgb)*globalparameter(itgbi)**2
2945  CALL addsum(adder)
2946  END DO
2947  END IF
2948 
2949 
2950  ! ----- add contributions from "measurement" -----------------------
2951 
2952 
2953  i=1
2954  DO WHILE (i <= lenmeasurements)
2955  rhs=listmeasurements(i )%value ! right hand side
2956  sgm=listmeasurements(i+1)%value ! sigma parameter
2957  i=i+2
2958  weight=0.0
2959  IF(sgm > 0.0) weight=1.0/sgm**2
2960 
2961  dsum=-rhs
2962 
2963  ! loop over label/factor pairs
2964  ia=i
2965  DO
2966  i=i+1
2967  IF(i > lenmeasurements) EXIT
2968  IF(listmeasurements(i)%label == 0) EXIT
2969  END DO
2970  ib=i-1
2971 
2972  DO j=ia,ib
2973  factrj=listmeasurements(j)%value
2974  itgbij=inone(listmeasurements(j)%label) ! total parameter index
2975  IF(itgbij /= 0) THEN
2976  dsum=dsum+factrj*globalparameter(itgbij) ! residuum
2977  END IF
2978  ! add to vector
2979  ivgbij=0
2980  IF(itgbij /= 0) ivgbij=globalparlabelindex(2,itgbij) ! variable-parameter index
2981  IF(ivgbij > 0) THEN
2982  globalvector(ivgbij)=globalvector(ivgbij) -weight*dsum*factrj ! vector
2983  globalcounter(ivgbij)=globalcounter(ivgbij)+1
2984  END IF
2985 
2986  IF(icalcm == 1.AND.ivgbij > 0) THEN
2987  DO k=ia,j
2988  factrk=listmeasurements(k)%value
2989  itgbik=inone(listmeasurements(k)%label) ! total parameter index
2990  ! add to matrix
2991  ivgbik=0
2992  IF(itgbik /= 0) ivgbik=globalparlabelindex(2,itgbik) ! variable-parameter index
2993  IF(ivgbij > 0.AND.ivgbik > 0) THEN !
2994  CALL mupdat(ivgbij,ivgbik,weight*factrj*factrk)
2995  END IF
2996  END DO
2997  END IF
2998  END DO
2999 
3000  adder=weight*dsum**2
3001  CALL addsum(adder) ! accumulate chi-square
3002 
3003  END DO
3004 
3005  ! ----- printout ---------------------------------------------------
3006 
3007 
3008  CALL getsum(fvalue) ! get accurate sum (Chi^2)
3009  flines=0.5_mpd*fvalue ! Likelihood function value
3010  rloop=iterat+0.01*nloopn
3011  actfun=real(funref-fvalue,mps)
3012  IF(nloopn == 1) actfun=0.0
3013  ngras=nint(angras,mpi)
3014  ratae=0.0 !!!
3015  IF(delfun /= 0.0) THEN
3016  ratae=min(99.9,actfun/delfun) !!!
3017  ratae=max(-99.9,ratae)
3018  END IF
3019 
3020  ! rejects ...
3021 
3022  nrej =nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3)
3023  IF(nloopn == 1) THEN
3024  IF(nrej /= 0) THEN
3025  WRITE(*,*) ' '
3026  WRITE(*,*) 'Data rejected in initial loop:'
3027  WRITE(*,*) ' ', &
3028  nrejec(0), ' (rank deficit/NaN) ',nrejec(1),' (Ndf=0) ', &
3029  nrejec(2), ' (huge) ',nrejec(3),' (large)'
3030  END IF
3031  END IF
3032  ! IF(NREJEC(1)+NREJEC(2)+NREJEC(3).NE.0) THEN
3033  ! WRITE(LUNLOG,*) 'Data rejected in initial loop:',NREJEC(1),
3034  ! + ' (Ndf=0) ',NREJEC(2),' (huge) ',NREJEC(3),' (large)'
3035  ! END IF
3036 
3037 
3038  IF(newite.AND.iterat == 2) THEN
3039  IF(nrecpr /= 0.OR.nrecp2 /= 0) nrecer=nrec3
3040  IF(nrecpr < 0) THEN
3041  nrecpr=nrec1
3042  END IF
3043  IF(nrecp2 < 0) THEN
3044  nrecp2=nrec2
3045  END IF
3046  END IF
3047 
3048  IF(nloopn <= 2) THEN
3049  IF(nhistp /= 0) THEN
3050  ! CALL HMPRNT(3) ! scaled residual of single measurement
3051  ! CALL HMPRNT(12) ! scaled residual of single measurement
3052  ! CALL HMPRNT(4) ! chi^2/Ndf
3053  END IF
3054  CALL hmpwrt(3)
3055  CALL hmpwrt(12)
3056  CALL hmpwrt(4)
3057  CALL gmpwrt(4) ! location, dispersion (res.) as a function of record nr
3058  IF (nloopn <= lfitnp) THEN
3059  CALL hmpwrt(13)
3060  CALL hmpwrt(14)
3061  CALL gmpwrt(5) ! location, dispersion (pull) as a function of record nr
3062  END IF
3063  END IF
3064  ! IF(NLOOPN.EQ.2.AND.NHISTP.NE.0) CALL HMPRNT(6)
3065  IF(nloopn == 2) CALL hmpwrt(6)
3066  IF(nloopn <= 1) THEN
3067  ! IF(NHISTP.NE.0) CALL HMPRNT(5) ! number of degrees of freedom
3068  ! IF(NHISTP.NE.0) CALL HMPRNT(11) ! Nlocal
3069  CALL hmpwrt(5)
3070  CALL hmpwrt(11)
3071  END IF
3072 
3073  ! local fit: band matrix structure !?
3074  IF (nloopn == 1.AND.nbndr(1)+nbndr(2) > 0) THEN
3075  DO lun=6,8,2
3076  WRITE(lun,*) ' '
3077  WRITE(lun,*) ' === local fits have bordered band matrix structure ==='
3078  IF (nbndr(1) > 0 ) WRITE(lun,101) ' NBNDR',nbndr(1),'number of records (upper/left border)'
3079  IF (nbndr(2) > 0 ) WRITE(lun,101) ' NBNDR',nbndr(2),'number of records (lower/right border)'
3080  WRITE(lun,101) ' NBDRX',nbdrx,'max border size'
3081  WRITE(lun,101) ' NBNDX',nbndx,'max band width'
3082  END DO
3083  END IF
3084 
3085  lastit=iterat
3086 
3087  ! monitoring of residuals
3088  IF (imonit < 0 .OR. (nloopn == 1 .AND. btest(imonit,0))) CALL monres
3089 
3090 101 FORMAT(1x,a8,' =',i10,' = ',a)
3091 ! 101 FORMAT(' LOOPN',I6,' Function value',F22.8,10X,I6,' records')
3092 ! 102 FORMAT(' incl. constraint penalty',F22.8)
3093 ! 103 FORMAT(I13,3X,A,G12.4)
3094 END SUBROUTINE loopn ! loop with fits
3095 
3099 
3100 SUBROUTINE ploopa(lunp)
3101  USE mpmod
3102 
3103  IMPLICIT NONE
3104 
3105  INTEGER(mpi), INTENT(IN) :: lunp
3106  ! ..
3107  WRITE(lunp,*) ' '
3108  WRITE(lunp,101) ! header line
3109  WRITE(lunp,102) ! header line
3110 101 FORMAT(' it fc',' fcn_value dfcn_exp slpr costh iit st', &
3111  ' ls step cutf',1x,'rejects hhmmss FMS')
3112 102 FORMAT(' -- --',' ----------- -------- ---- ----- --- --', &
3113  ' -- ----- ----',1x,'------- ------ ---')
3114  RETURN
3115 END SUBROUTINE ploopa ! title for iteration
3116 
3120 
3121 SUBROUTINE ploopb(lunp)
3122  USE mpmod
3123 
3124  IMPLICIT NONE
3125  INTEGER(mpi) :: ma
3126  INTEGER :: minut
3127  INTEGER(mpi) :: nfa
3128  INTEGER :: nhour
3129  INTEGER(mpi) :: nrej
3130  INTEGER(mpi) :: nsecnd
3131  REAL(mps) :: ratae
3132  REAL :: rstb
3133  REAL(mps) :: secnd
3134  REAL(mps) :: slopes(3)
3135  REAL(mps) :: steps(3)
3136  REAL, DIMENSION(2) :: ta
3137  REAl etime
3138 
3139  INTEGER(mpi), INTENT(IN) :: lunp
3140 
3141  CHARACTER (LEN=4):: ccalcm(4)
3142  DATA ccalcm / ' end',' S', ' F ',' FMS' /
3143  SAVE
3144 
3145  nrej=nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) ! rejects
3146  IF(nrej > 9999999) nrej=9999999
3147  rstb=etime(ta)
3148  deltim=rstb-rstart
3149  CALL sechms(deltim,nhour,minut,secnd) ! time
3150  nsecnd=nint(secnd,mpi)
3151  IF(iterat == 0) THEN
3152  WRITE(lunp,103) iterat,nloopn,fvalue, &
3153  chicut,nrej,nhour,minut,nsecnd,ccalcm(lcalcm)
3154  ELSE
3155  IF (lsinfo == 10) THEN ! line search skipped
3156  WRITE(lunp,105) iterat,nloopn,fvalue,delfun, &
3157  iitera,istopa,chicut,nrej,nhour,minut,nsecnd,ccalcm(lcalcm)
3158  ELSE
3159  CALL ptlopt(nfa,ma,slopes,steps) ! slopes steps
3160  ratae=max(-99.9,min(99.9,slopes(2)/slopes(1)))
3161  stepl=steps(2)
3162  WRITE(lunp,104) iterat,nloopn,fvalue,delfun,ratae,angras, &
3163  iitera,istopa,lsinfo,stepl, chicut,nrej,nhour,minut,nsecnd,ccalcm(lcalcm)
3164  ENDIF
3165  END IF
3166 103 FORMAT(i3,i3,e12.5,38x,f5.1, 1x,i7, i3,i2.2,i2.2,a4)
3167 104 FORMAT(i3,i3,e12.5,1x,e8.2,f6.3,f6.3,i5,2i3,f6.3,f5.1, &
3168  1x,i7, i3,i2.2,i2.2,a4)
3169 105 FORMAT(i3,i3,e12.5,1x,e8.2,12x,i5,i3,9x,f5.1, &
3170  1x,i7, i3,i2.2,i2.2,a4)
3171  RETURN
3172 END SUBROUTINE ploopb ! iteration line
3173 
3177 
3178 SUBROUTINE ploopc(lunp)
3179  USE mpmod
3180 
3181  IMPLICIT NONE
3182  INTEGER(mpi) :: ma
3183  INTEGER(mpi) :: minut
3184  INTEGER(mpi) :: nfa
3185  INTEGER(mpi) :: nhour
3186  INTEGER(mpi) :: nrej
3187  INTEGER(mpi) :: nsecnd
3188  REAL(mps) :: ratae
3189  REAL :: rstb
3190  REAL(mps) :: secnd
3191  REAL(mps) :: slopes(3)
3192  REAL(mps) :: steps(3)
3193  REAL, DIMENSION(2) :: ta
3194  REAL etime
3195 
3196  INTEGER(mpi), INTENT(IN) :: lunp
3197  CHARACTER (LEN=4):: ccalcm(4)
3198  DATA ccalcm / ' end',' S', ' F ',' FMS' /
3199  SAVE
3200 
3201  nrej=nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) ! rejects
3202  IF(nrej > 9999999) nrej=9999999
3203  rstb=etime(ta)
3204  deltim=rstb-rstart
3205  CALL sechms(deltim,nhour,minut,secnd) ! time
3206  nsecnd=nint(secnd,mpi)
3207  IF (lsinfo == 10) THEN ! line search skipped
3208  WRITE(lunp,104) nloopn,fvalue,nrej,nhour,minut,nsecnd,ccalcm(lcalcm)
3209  ELSE
3210  CALL ptlopt(nfa,ma,slopes,steps) ! slopes steps
3211  ratae=abs(slopes(2)/slopes(1))
3212  stepl=steps(2)
3213  WRITE(lunp,105) nloopn,fvalue, ratae,lsinfo, &
3214  stepl,nrej,nhour,minut,nsecnd,ccalcm(lcalcm)
3215  END IF
3216 104 FORMAT(3x,i3,e12.5,9x, 35x, i7, i3,i2.2,i2.2,a4)
3217 105 FORMAT(3x,i3,e12.5,9x, f6.3,14x,i3,f6.3,6x, i7, i3,i2.2,i2.2,a4)
3218  RETURN
3219 
3220 END SUBROUTINE ploopc ! sub-iteration line
3221 
3225 
3226 SUBROUTINE ploopd(lunp)
3227  USE mpmod
3228  IMPLICIT NONE
3229  INTEGER :: minut
3230  INTEGER :: nhour
3231  INTEGER(mpi) :: nsecnd
3232  REAL :: rstb
3233  REAL(mps) :: secnd
3234  REAL, DIMENSION(2) :: ta
3235  REAL etime
3236 
3237  INTEGER(mpi), INTENT(IN) :: lunp
3238  CHARACTER (LEN=4):: ccalcm(4)
3239  DATA ccalcm / ' end',' S', ' F ',' FMS' /
3240  SAVE
3241  rstb=etime(ta)
3242  deltim=rstb-rstart
3243  CALL sechms(deltim,nhour,minut,secnd) ! time
3244  nsecnd=nint(secnd,mpi)
3245 
3246  WRITE(lunp,106) nhour,minut,nsecnd,ccalcm(lcalcm)
3247 106 FORMAT(69x,i3,i2.2,i2.2,a4)
3248  RETURN
3249 END SUBROUTINE ploopd
3250 
3252 SUBROUTINE explfc(lunit)
3253  USE mpdef
3254  USE mpmod, ONLY: metsol
3255 
3256  IMPLICIT NONE
3257  INTEGER(mpi) :: lunit
3258  WRITE(lunit,*) ' '
3259  WRITE(lunit,102) 'Explanation of iteration table'
3260  WRITE(lunit,102) '=============================='
3261  WRITE(lunit,101) 'it', &
3262  'iteration number. Global parameters are improved for it > 0.'
3263  WRITE(lunit,102) 'First function evaluation is called iteraton 0.'
3264  WRITE(lunit,101) 'fc', 'number of function evaluations.'
3265  WRITE(lunit,101) 'fcn_value', 'value of 2 x Likelihood function (LF).'
3266  WRITE(lunit,102) 'The final value is the chi^2 value of the fit and should'
3267  WRITE(lunit,102) 'be about equal to the NDF (see below).'
3268  WRITE(lunit,101) 'dfcn_exp', &
3269  'expected reduction of the value of the Likelihood function (LF)'
3270  WRITE(lunit,101) 'slpr', 'ratio of the actual slope to inital slope.'
3271  WRITE(lunit,101) 'costh', &
3272  'cosine of angle between search direction and -gradient'
3273  IF (metsol == 4) THEN
3274  WRITE(lunit,101) 'iit', &
3275  'number of internal iterations in MINRES algorithm'
3276  WRITE(lunit,101) 'st', 'stop code of MINRES algorithm'
3277  WRITE(lunit,102) '< 0: rhs is very special, with beta2 = 0'
3278  WRITE(lunit,102) '= 0: rhs b = 0, i.e. the exact solution is x = 0'
3279  WRITE(lunit,102) '= 1 requested accuracy achieved, as determined by rtol'
3280  WRITE(lunit,102) '= 2 reasonable accuracy achieved, given eps'
3281  WRITE(lunit,102) '= 3 x has converged to an eigenvector'
3282  WRITE(lunit,102) '= 4 matrix ill-conditioned (Acond has exceeded 0.1/eps)'
3283  WRITE(lunit,102) '= 5 the iteration limit was reached'
3284  WRITE(lunit,102) '= 6 Matrix x vector does not define a symmetric matrix'
3285  WRITE(lunit,102) '= 7 Preconditioner does not define a symmetric matrix'
3286  ELSEIF (metsol == 5) THEN
3287  WRITE(lunit,101) 'iit', &
3288  'number of internal iterations in MINRES-QLP algorithm'
3289  WRITE(lunit,101) 'st', 'stop code of MINRES-QLP algorithm'
3290  WRITE(lunit,102) '= 1: beta_{k+1} < eps, iteration k is the final Lanczos step.'
3291  WRITE(lunit,102) '= 2: beta2 = 0. If M = I, b and x are eigenvectors of A.'
3292  WRITE(lunit,102) '= 3: beta1 = 0. The exact solution is x = 0.'
3293  WRITE(lunit,102) '= 4: A solution to (poss. singular) Ax = b found, given rtol.'
3294  WRITE(lunit,102) '= 5: A solution to (poss. singular) Ax = b found, given eps.'
3295  WRITE(lunit,102) '= 6: Pseudoinverse solution for singular LS problem, given rtol.'
3296  WRITE(lunit,102) '= 7: Pseudoinverse solution for singular LS problem, given eps.'
3297  WRITE(lunit,102) '= 8: The iteration limit was reached.'
3298  WRITE(lunit,102) '= 9: The operator defined by Aprod appears to be unsymmetric.'
3299  WRITE(lunit,102) '=10: The operator defined by Msolve appears to be unsymmetric.'
3300  WRITE(lunit,102) '=11: The operator defined by Msolve appears to be indefinite.'
3301  WRITE(lunit,102) '=12: xnorm has exceeded maxxnorm or will exceed it next iteration.'
3302  WRITE(lunit,102) '=13: Acond has exceeded Acondlim or 0.1/eps.'
3303  WRITE(lunit,102) '=14: Least-squares problem but no converged solution yet.'
3304  WRITE(lunit,102) '=15: A null vector obtained, given rtol.'
3305  ENDIF
3306  WRITE(lunit,101) 'ls', 'line search info'
3307  WRITE(lunit,102) '< 0 recalculate function'
3308  WRITE(lunit,102) '= 0: N or STP lt 0 or step not descending'
3309  WRITE(lunit,102) '= 1: Linesearch convergence conditions reached'
3310  WRITE(lunit,102) '= 2: interval of uncertainty at lower limit'
3311  WRITE(lunit,102) '= 3: max nr of line search calls reached'
3312  WRITE(lunit,102) '= 4: step at the lower bound'
3313  WRITE(lunit,102) '= 5: step at the upper bound'
3314  WRITE(lunit,102) '= 6: rounding error limitation'
3315  WRITE(lunit,101) 'step', &
3316  'the factor for the Newton step during the line search. Usually'
3317  WRITE(lunit,102) &
3318  'a value of 1 gives a sufficient reduction of the LF. Oherwise'
3319  WRITE(lunit,102) 'other step values are tried.'
3320  WRITE(lunit,101) 'cutf', &
3321  'cut factor. Local fits are rejected, if their chi^2 value'
3322  WRITE(lunit,102) &
3323  'is larger than the 3-sigma chi^2 value times the cut factor.'
3324  WRITE(lunit,102) 'A cut factor of 1 is used finally, but initially a larger'
3325  WRITE(lunit,102) 'factor may be used. A value of 0.0 means no cut.'
3326  WRITE(lunit,101) 'rejects', 'total number of rejected local fits.'
3327  WRITE(lunit,101) 'hmmsec', 'the time in hours (h), minutes (mm) and seconds.'
3328  WRITE(lunit,101) 'FMS', 'calculation of Function value, Matrix, Solution.'
3329  WRITE(lunit,*) ' '
3330 
3331 101 FORMAT(a9,' = ',a)
3332 102 FORMAT(13x,a)
3333 END SUBROUTINE explfc
3334 
3342 
3343 SUBROUTINE mupdat(i,j,add) !
3344  USE mpmod
3345 
3346  IMPLICIT NONE
3347 
3348  INTEGER(mpi), INTENT(IN) :: i
3349  INTEGER(mpi), INTENT(IN) :: j
3350  REAL(mpd), INTENT(IN) :: add
3351 
3352  INTEGER(mpl):: ijadd
3353  INTEGER(mpl):: ia
3354  INTEGER(mpl):: ja
3355  INTEGER(mpl):: ij
3356  ! ...
3357  IF(i <= 0.OR.j <= 0) RETURN
3358  ia=max(i,j) ! larger
3359  ja=min(i,j) ! smaller
3360  ij=0
3361  IF(matsto == 2) THEN ! sparse symmetric matrix
3362  ij=ijadd(i,j) ! inline code requires same time
3363  IF (ij == 0) RETURN ! pair is suppressed
3364  IF (ij > 0) THEN
3365  globalmatd(ij)=globalmatd(ij)+add
3366  ELSE
3367  globalmatf(-ij)=globalmatf(-ij)+real(add,mps)
3368  END IF
3369  ELSE ! full or unpacked (block diagonal) symmetric matrix
3370  ! global (ia,ib) to local (row,col) in block
3371  ij=globalrowoffsets(ia)+ja
3372  globalmatd(ij)=globalmatd(ij)+add
3373  END IF
3374  IF(metsol >= 4.AND.metsol < 7) THEN
3375  IF(mbandw > 0) THEN ! for Cholesky decomposition
3376  IF(ia <= nvgb) THEN ! variable global parameter
3377  ij=indprecond(ia)-ia+ja
3378  IF(ia > 1.AND.ij <= indprecond(ia-1)) ij=0
3379  IF(ij /= 0) matprecond(ij)=matprecond(ij)+add
3380  IF(ij < 0.OR.ij > size(matprecond)) THEN
3381  CALL peend(23,'Aborted, bad matrix index')
3382  stop 'mupdat: bad index'
3383  END IF
3384  ELSE ! Lagrange multiplier
3385  ij=indprecond(nvgb)+(ia-nvgb-1)*nvgb+ja
3386  IF(ij /= 0) matprecond(ij)=matprecond(ij)+add
3387  END IF
3388  ELSE IF(mbandw == 0) THEN ! default preconditioner
3389  IF(ia <= nvgb) THEN ! variable global parameter
3390  IF(ja == ia) matprecond(ia)=matprecond(ia)+add ! diag
3391  ELSE ! Lagrange multiplier
3392  ij=nvgb+(ia-nvgb-1)*nvgb+ja
3393  IF(ij /= 0) matprecond(ij)=matprecond(ij)+add
3394  END IF
3395  END IF
3396  END IF
3397 END SUBROUTINE mupdat
3398 
3399 
3410 
3411 SUBROUTINE mgupdt(i,j1,j2,il,jl,sub)
3412  USE mpmod
3413 
3414  IMPLICIT NONE
3415 
3416  INTEGER(mpi), INTENT(IN) :: i
3417  INTEGER(mpi), INTENT(IN) :: j1
3418  INTEGER(mpi), INTENT(IN) :: j2
3419  INTEGER(mpi), INTENT(IN) :: il
3420  INTEGER(mpi), INTENT(IN) :: jl
3421  REAL(mpd), INTENT(IN) :: sub(*)
3422 
3423  INTEGER(mpl):: ij
3424  INTEGER(mpi):: ia
3425  INTEGER(mpi):: ib
3426  INTEGER(mpi):: iprc
3427  INTEGER(mpi):: ir
3428  INTEGER(mpi):: ja
3429  INTEGER(mpi):: jb
3430  INTEGER(mpi):: ijl
3431  INTEGER(mpi):: k
3432  INTEGER(mpi):: lr
3433  INTEGER(mpi):: nc
3434  ! ...
3435  IF(i <= 0.OR.j1 <= 0.OR.j2 > i) RETURN
3436  ia=globalallindexgroups(i) ! first (global) row
3437  ib=globalallindexgroups(i+1)-1 ! last (global) row
3438  ja=globalallindexgroups(j1) ! first (global) column
3439  jb=globalallindexgroups(j2+1)-1 ! last (global) column
3440 
3441  IF(matsto == 2) THEN ! sparse symmetric matrix
3442  CALL ijpgrp(i,j1,ij,lr,iprc) ! index of first element of group 'j1'
3443  !print *, ' mgupdt ', i,j1,j2,il,jl,ij,lr,iprc
3444  IF (ij == 0) THEN
3445  print *, ' MGUPDT: ij=0', i,j1,j2,il,jl,ij,lr,iprc
3446  stop
3447  END IF
3448  k=il
3449  ijl=(k*k-k)/2 ! ISYM index offset (subtrahends matrix)
3450  DO ir=ia,ib
3451  nc=min(ir,jb)-ja ! number of columns -1
3452  IF (jb >= ir) THEN ! diagonal element
3453  globalmatd(ir)=globalmatd(ir)-sub(ijl+jl+nc)
3454  nc=nc-1
3455  END IF
3456  ! off-diagonal elements
3457  IF (iprc == 1) THEN
3458  globalmatd(ij:ij+nc)=globalmatd(ij:ij+nc)-sub(ijl+jl:ijl+jl+nc)
3459  ELSE
3460  globalmatf(ij:ij+nc)=globalmatf(ij:ij+nc)-real(sub(ijl+jl:ijl+jl+nc),mps)
3461  END IF
3462  ij=ij+lr
3463  ijl=ijl+k
3464  k=k+1
3465  END DO
3466  ELSE ! full or unpacked (block diagonal) symmetric matrix
3467  k=il
3468  ijl=(k*k-k)/2 ! ISYM index offset (subtrahends matrix)
3469  DO ir=ia,ib
3470  ! global (ir,0) to local (row,col) in block
3471  ij=globalrowoffsets(ir)
3472  nc=min(ir,jb)-ja ! number of columns -1
3473  globalmatd(ij+ja:ij+ja+nc)=globalmatd(ij+ja:ij+ja+nc)-sub(ijl+jl:ijl+jl+nc)
3474  ijl=ijl+k
3475  k=k+1
3476  END DO
3477  END IF
3478 
3479 END SUBROUTINE mgupdt
3480 
3481 
3511 
3512 SUBROUTINE loopbf(nrej,ndfs,sndf,dchi2s, numfil,naccf,chi2f,ndff)
3513  USE mpmod
3514 
3515  IMPLICIT NONE
3516  REAL(mpd) :: cauchy
3517  REAL(mps) :: chichi
3518  REAL(mps) :: chlimt
3519  REAL(mps) :: chndf
3520  REAL(mpd) :: chuber
3521  REAL(mpd) :: down
3522  REAL(mpr8) :: glder
3523  REAL(mpd) :: pull
3524  REAL(mpd) :: r1
3525  REAL(mpd) :: r2
3526  REAL(mps) :: rec
3527  REAL(mpd) :: rerr
3528  REAL(mpd) :: resid
3529  REAL(mps) :: resing
3530  REAL(mpd) :: resmax
3531  REAL(mpd) :: rmeas
3532  REAL(mpd) :: rmloc
3533  REAL(mpd) :: suwt
3534  REAL(mps) :: used
3535  REAL(mpd) :: wght
3536  REAL(mps) :: chindl
3537  INTEGER(mpi) :: i
3538  INTEGER(mpi) :: ia
3539  INTEGER(mpi) :: ib
3540  INTEGER(mpi) :: ibuf
3541  INTEGER(mpi) :: ichunk
3542  INTEGER(mpl) :: icmn
3543  INTEGER(mpl) :: icost
3544  INTEGER(mpi) :: id
3545  INTEGER(mpi) :: idiag
3546  INTEGER(mpi) :: ieq
3547  INTEGER(mpi) :: iext
3548  INTEGER(mpi) :: ij
3549  INTEGER(mpi) :: ije
3550  INTEGER(mpi) :: ijn
3551  INTEGER(mpi) :: ijsym
3552  INTEGER(mpi) :: ik
3553  INTEGER(mpi) :: ike
3554  INTEGER(mpi) :: il
3555  INTEGER(mpi) :: im
3556  INTEGER(mpi) :: imeas
3557  INTEGER(mpi) :: in
3558  INTEGER(mpi) :: inder
3559  INTEGER(mpi) :: inv
3560  INTEGER(mpi) :: ioffb
3561  INTEGER(mpi) :: ioffc
3562  INTEGER(mpi) :: ioffd
3563  INTEGER(mpi) :: ioffe
3564  INTEGER(mpi) :: ioffi
3565  INTEGER(mpi) :: ioffq
3566  INTEGER(mpi) :: iprc
3567  INTEGER(mpi) :: iprcnx
3568  INTEGER(mpi) :: iprdbg
3569  INTEGER(mpi) :: iproc
3570  INTEGER(mpi) :: irbin
3571  INTEGER(mpi) :: irow
3572  INTEGER(mpi) :: isfrst
3573  INTEGER(mpi) :: isize
3574  INTEGER(mpi) :: islast
3575  INTEGER(mpi) :: ist
3576  INTEGER(mpi) :: iter
3577  INTEGER(mpi) :: itgbi
3578  INTEGER(mpi) :: ivgbj
3579  INTEGER(mpi) :: ivgbk
3580  INTEGER(mpi) :: ivpgrp
3581  INTEGER(mpi) :: j
3582  INTEGER(mpi) :: j1
3583  INTEGER(mpi) :: ja
3584  INTEGER(mpi) :: jb
3585  INTEGER(mpi) :: jk
3586  INTEGER(mpi) :: jl
3587  INTEGER(mpi) :: jl1
3588  INTEGER(mpi) :: jn
3589  INTEGER(mpi) :: jnx
3590  INTEGER(mpi) :: joffd
3591  INTEGER(mpi) :: joffi
3592  INTEGER(mpi) :: jproc
3593  INTEGER(mpi) :: jsp
3594  INTEGER(mpi) :: k
3595  INTEGER(mpi) :: kbdr
3596  INTEGER(mpi) :: kbdrx
3597  INTEGER(mpi) :: kbnd
3598  INTEGER(mpi) :: kfl
3599  INTEGER(mpi) :: kx
3600  INTEGER(mpi) :: lvpgrp
3601  INTEGER(mpi) :: mbdr
3602  INTEGER(mpi) :: mbnd
3603  INTEGER(mpi) :: mside
3604  INTEGER(mpi) :: nalc
3605  INTEGER(mpi) :: nalg
3606  INTEGER(mpi) :: nan
3607  INTEGER(mpi) :: nb
3608  INTEGER(mpi) :: ndf
3609  INTEGER(mpi) :: ndfsd
3610  INTEGER(mpi) :: ndown
3611  INTEGER(mpi) :: neq
3612  INTEGER(mpi) :: nfred
3613  INTEGER(mpi) :: nfrei
3614  INTEGER(mpi) :: ngg
3615  INTEGER(mpi) :: nprdbg
3616  INTEGER(mpi) :: nrank
3617  INTEGER(mpi) :: nrc
3618  INTEGER(mpi) :: nst
3619  INTEGER(mpi) :: nter
3620  INTEGER(mpi) :: nweig
3621  INTEGER(mpi) :: ngrp
3622  INTEGER(mpi) :: npar
3623 
3624  INTEGER(mpi), INTENT(IN OUT) :: nrej(0:3)
3625  INTEGER(mpi), INTENT(IN OUT) :: ndfs
3626  REAL(mpd), INTENT(IN OUT) :: sndf
3627  REAL(mpd), INTENT(IN OUT) :: dchi2s
3628  INTEGER(mpi), INTENT(IN) :: numfil
3629  INTEGER(mpi), INTENT(IN OUT) :: naccf(numfil)
3630  REAL(mps), INTENT(IN OUT) :: chi2f(numfil)
3631  INTEGER(mpi), INTENT(IN OUT) :: ndff(numfil)
3632 
3633  REAL(mpd):: dchi2
3634  REAL(mpd)::dvar
3635  REAL(mpd):: dw1
3636  REAL(mpd)::dw2
3637  REAL(mpd)::summ
3638  INTEGER(mpi) :: ijprec
3639 
3640  !$ INTEGER(mpi) OMP_GET_THREAD_NUM
3641 
3642  LOGICAL:: lprnt
3643  LOGICAL::lhist
3644 
3645  CHARACTER (LEN=3):: chast
3646  DATA chuber/1.345_mpd/ ! constant for Huber down-weighting
3647  DATA cauchy/2.3849_mpd/ ! constant for Cauchy down-weighting
3648  SAVE chuber,cauchy
3649  ! ...
3650  ijsym(i,j)=min(i,j)+(max(i,j)*max(i,j)-max(i,j))/2
3651  isfrst(ibuf)=readbufferpointer(ibuf)+1
3652  islast(ibuf)=readbufferdatai(readbufferpointer(ibuf))
3653  inder(i)=readbufferdatai(i)
3654  glder(i)=readbufferdatad(i)
3655 
3656  ichunk=min((numreadbuffer+mthrd-1)/mthrd/32+1,256)
3657  ! reset header, 3 words per thread:
3658  ! number of entries, offset to data, indices
3659  writebufferinfo=0
3660  writebufferdata=0.
3661  nprdbg=0
3662  iprdbg=-1
3663 
3664  ! parallelize record loop
3665  ! private copy of NDFS,.. for each thread, combined at end, init with 0.
3666  !$OMP PARALLEL DO &
3667  !$OMP DEFAULT(PRIVATE) &
3668  !$OMP SHARED(numReadBuffer,readBufferPointer,readBufferDataI, &
3669  !$OMP readBufferDataD,writeBufferHeader,writeBufferInfo, &
3670  !$OMP writeBufferData,writeBufferIndices,writeBufferUpdates,globalVector,globalCounter, &
3671  !$OMP globalParameter,globalParLabelIndex,globalIndexUsage,backIndexUsage, &
3672  !$OMP measBins,numMeas,measIndex,measRes,measHists,globalAllParToGroup,globalAllIndexGroups, &
3673  !$OMP localCorrections,localEquations, &
3674  !$OMP NAGB,NVGB,NAGBN,ICALCM,ICHUNK,NLOOPN,NRECER,NPRDBG,IPRDBG, &
3675  !$OMP NEWITE,CHICUT,LHUBER,CHUBER,ITERAT,NRECPR,MTHRD,NSPC,NAEQN, &
3676  !$OMP DWCUT,CHHUGE,NRECP2,CAUCHY,LFITNP,LFITBB,IMONIT,IMONMD,MONPG1,LUNLOG) &
3677  !$OMP REDUCTION(+:NDFS,SNDF,DCHI2S,NREJ,NBNDR,NACCF,CHI2F,NDFF) &
3678  !$OMP REDUCTION(MAX:NBNDX,NBDRX) &
3679  !$OMP REDUCTION(MIN:NREC3) &
3680  !$OMP SCHEDULE(DYNAMIC,ICHUNK)
3681  DO ibuf=1,numreadbuffer ! buffer for current record
3682  nrc=readbufferdatai(isfrst(ibuf)-2) ! record
3683  kfl=nint(readbufferdatad(isfrst(ibuf)-1),mpi) ! file
3684  dw1=real(readbufferdatad(isfrst(ibuf)-2),mpd) ! weight
3685  dw2=sqrt(dw1)
3686 
3687  iproc=0
3688  !$ IPROC=OMP_GET_THREAD_NUM() ! thread number
3689  ioffb=nagb*iproc ! offset 'f'.
3690  ioffc=nagbn*iproc ! offset 'c'.
3691  ioffe=nvgb*iproc ! offset 'e'
3692  ioffd=writebufferheader(-1)*iproc+writebufferinfo(2,iproc+1) ! offset data
3693  ioffi=writebufferheader(1)*iproc+writebufferinfo(3,iproc+1)+2 ! offset indices
3694  ioffq=naeqn*iproc ! offset equations (measurements)
3695  ! ----- reset ------------------------------------------------------
3696  lprnt=.false.
3697  lhist=(iproc == 0)
3698  rec=nrc ! floating point value
3699  IF(nloopn == 1.AND.mod(nrc,100000) == 0) THEN
3700  WRITE(*,*) 'Record',nrc,' ... still reading'
3701  IF(monpg1>0) WRITE(lunlog,*) 'Record',nrc,' ... still reading'
3702  END IF
3703 
3704  ! printout/debug only for one thread at a time
3705 
3706 
3707  ! flag for record printout -----------------------------------------
3708 
3709  lprnt=.false.
3710  IF(newite.AND.(iterat == 1.OR.iterat == 3)) THEN
3711  IF(nrc == nrecpr) lprnt=.true.
3712  IF(nrc == nrecp2) lprnt=.true.
3713  IF(nrc == nrecer) lprnt=.true.
3714  END IF
3715  IF (lprnt)THEN
3716  !$OMP ATOMIC
3717  nprdbg=nprdbg+1 ! number of threads with debug
3718  IF (nprdbg == 1) iprdbg=iproc ! first thread with debug
3719  IF (iproc /= iprdbg) lprnt=.false.
3720  ! print *, ' LPRNT ', NRC, NPRDBG, IPRDBG, IPROC, LPRNT
3721  END IF
3722  IF(lprnt) THEN
3723  WRITE(1,*) ' '
3724  WRITE(1,*) '------------------ Loop',nloopn, &
3725  ': Printout for record',nrc,iproc
3726  WRITE(1,*) ' '
3727  END IF
3728 
3729  ! ----- print data -------------------------------------------------
3730 
3731  IF(lprnt) THEN
3732  imeas=0 ! local derivatives
3733  ist=isfrst(ibuf)
3734  nst=islast(ibuf)
3735  DO ! loop over measurements
3736  CALL isjajb(nst,ist,ja,jb,jsp)
3737  IF(ja == 0) EXIT
3738  IF(imeas == 0) WRITE(1,1121)
3739  imeas=imeas+1
3740  WRITE(1,1122) imeas,glder(ja),glder(jb), &
3741  (inder(ja+j),glder(ja+j),j=1,jb-ja-1)
3742  END DO
3743 1121 FORMAT(/'Measured value and local derivatives'/ &
3744  ' i measured std_dev index...derivative ...')
3745 1122 FORMAT(i3,2g12.4,3(i3,g12.4)/(27x,3(i3,g12.4)))
3746 
3747  imeas=0 ! global derivatives
3748  ist=isfrst(ibuf)
3749  nst=islast(ibuf)
3750  DO ! loop over measurements
3751  CALL isjajb(nst,ist,ja,jb,jsp)
3752  IF(ja == 0) EXIT
3753  IF(imeas == 0) WRITE(1,1123)
3754  imeas=imeas+1
3755  IF (jb < ist) THEN
3756  IF(ist-jb > 2) THEN
3757  WRITE(1,1124) imeas,(globalparlabelindex(1,inder(jb+j)),inder(jb+j), &
3758  globalparlabelindex(2,inder(jb+j)),glder(jb+j),j=1,ist-jb)
3759  ELSE
3760  WRITE(1,1125) imeas,(globalparlabelindex(1,inder(jb+j)),inder(jb+j), &
3761  globalparlabelindex(2,inder(jb+j)),glder(jb+j),j=1,ist-jb)
3762  END IF
3763  END IF
3764  END DO
3765 1123 FORMAT(/'Global derivatives'/ &
3766  ' i label gindex vindex derivative ...')
3767 1124 FORMAT(i3,2(i9,i7,i7,g12.4)/(3x,2(i9,i7,i7,g12.4)))
3768 1125 FORMAT(i3,2(i9,i7,i7,g12.4))
3769  END IF
3770 
3771  ! ----- first loop -------------------------------------------------
3772  ! ------ prepare local fit ------
3773  ! count local and global derivates
3774  ! subtract actual alignment parameters from the measured data
3775 
3776  IF(lprnt) THEN
3777  WRITE(1,*) ' '
3778  WRITE(1,*) 'Data corrections using values of global parameters'
3779  WRITE(1,*) '=================================================='
3780  WRITE(1,101)
3781  END IF
3782  nalg=0 ! count number of global derivatives
3783  nalc=0 ! count number of local derivatives
3784  neq=0 ! count number of equations
3785 
3786  ist=isfrst(ibuf)
3787  nst=islast(ibuf)
3788  DO ! loop over measurements
3789  CALL isjajb(nst,ist,ja,jb,jsp)
3790  IF(ja == 0) EXIT
3791  rmeas=real(glder(ja),mpd) ! data
3792  neq=neq+1 ! count equation
3793  localequations(1,ioffq+neq)=ja
3794  localequations(2,ioffq+neq)=jb
3795  localequations(3,ioffq+neq)=ist
3796  ! subtract global ... from measured value
3797  DO j=1,ist-jb ! global parameter loop
3798  itgbi=inder(jb+j) ! global parameter label
3799  rmeas=rmeas-real(glder(jb+j),mpd)*globalparameter(itgbi) ! subtract !!! reversed
3800  IF (icalcm == 1) THEN
3801  ij=globalparlabelindex(2,itgbi) ! index of variable global parameter
3802  IF(ij > 0) THEN
3803  ijn=backindexusage(ioffe+ij) ! get index of index
3804  IF(ijn == 0) THEN ! not yet included
3805  nalg=nalg+1 ! count
3806  globalindexusage(ioffc+nalg)=ij ! store global index
3807  backindexusage(ioffe+ij)=nalg ! store back index
3808  END IF
3809  END IF
3810  END IF
3811  END DO
3812  IF(lprnt) THEN
3813  IF (jb < ist) WRITE(1,102) neq,glder(ja),rmeas,glder(jb)
3814  END IF
3815  readbufferdatad(ja)=real(rmeas,mpr8) ! global contribution subtracted
3816  DO j=1,jb-ja-1 ! local parameter loop
3817  ij=inder(ja+j)
3818  nalc=max(nalc,ij) ! number of local parameters
3819  END DO
3820  END DO
3821 101 FORMAT(' index measvalue corrvalue sigma')
3822 102 FORMAT(i6,2x,2g12.4,' +-',g12.4)
3823 
3824  IF(nalc <= 0) GO TO 90
3825 
3826  ngg=(nalg*nalg+nalg)/2
3827  ngrp=0
3828  IF (icalcm == 1) THEN
3829  localglobalmatrix(:nalg*nalc)=0.0_mpd ! reset global-local matrix
3830  localglobalmap(:nalg*nalc)=0 ! reset global-local map
3831  ! store parameter group indices
3832  CALL sort1k(globalindexusage(ioffc+1),nalg) ! sort global par.
3833  lvpgrp=-1
3834  npar=0
3835  DO k=1,nalg
3836  iext=globalindexusage(ioffc+k)
3837  backindexusage(ioffe+iext)=k ! update back index
3838  ivpgrp=globalallpartogroup(iext) ! group
3839  IF (ivpgrp /= lvpgrp) THEN
3840  ngrp=ngrp+1
3841  writebufferindices(ioffi+ngrp)=ivpgrp ! global par group indices
3842  lvpgrp=ivpgrp
3843  npar=npar+globalallindexgroups(ivpgrp+1)-globalallindexgroups(ivpgrp)
3844  END IF
3845  END DO
3846  ! check NPAR==NALG
3847  IF (npar /= nalg) THEN
3848  print *, ' mismatch of number of global parameters ', nrc, nalg, npar, ngrp
3849  print *, globalindexusage(ioffc+1:ioffc+nalg)
3850  print *, writebufferindices(ioffi+1:ioffi+ngrp)
3851  j=0
3852  DO k=1,ngrp
3853  ivpgrp=writebufferindices(ioffi+k)
3854  j=j+globalallindexgroups(ivpgrp+1)-globalallindexgroups(ivpgrp)
3855  IF (globalallpartogroup(globalindexusage(ioffc+j)) /= ivpgrp) &
3856  print *, ' bad group ', k, j, ivpgrp, globalindexusage(ioffc+j)
3857  END DO
3858  stop ' mismatch of number of global parameters '
3859  ENDIF
3860  writebufferindices(ioffi-1)=nrc ! index header:
3861  writebufferindices(ioffi )=ngrp ! event number, number of global par groups
3862  DO k=1,ngg
3863  writebufferupdates(ioffd+k)=0.0_mpd ! reset global-global matrix
3864  END DO
3865  END IF
3866  ! ----- iteration start and check ---------------------------------
3867 
3868  nter=1 ! first loop without down-weighting
3869  IF(nloopn /= 1.AND.lhuber /= 0) nter=lhuber
3870  localcorrections(ioffq+1:ioffq+neq) = 0._mpd
3871 
3872  ! check matrix for bordered band structure (MBDR+MBND+1 <= NALC)
3873  mbnd=-1
3874  mbdr=nalc
3875  mside=-1 ! side (1: upper/left border, 2: lower/right border)
3876  DO i=1, 2*nalc
3877  ibandh(i)=0
3878  END DO
3879  irow=1
3880  idiag=1
3881  ndfsd=0
3882 
3883  iter=0
3884  resmax=0.0
3885  DO WHILE(iter < nter) ! outlier suppresssion iteration loop
3886  iter=iter+1
3887  resmax=0.0
3888  IF(lprnt) THEN
3889  WRITE(1,*) ' '
3890  WRITE(1,*) 'Outlier-suppression iteration',iter,' of',nter
3891  WRITE(1,*) '=========================================='
3892  WRITE(1,*) ' '
3893  imeas=0
3894  END IF
3895 
3896  ! ----- second loop ------------------------------------------------
3897  ! accumulate normal equations for local fit and determine solution
3898  DO i=1,nalc
3899  blvec(i)=0.0_mpd ! reset vector
3900  END DO
3901  DO i=1,(nalc*nalc+nalc)/2 ! GF: FIXME - not really, local parameter number...
3902  clmat(i)=0.0_mpd ! (p)reset matrix
3903  END DO
3904  ndown=0
3905  nweig=0
3906  DO ieq=1,neq! loop over measurements
3907  ja=localequations(1,ioffq+ieq)
3908  jb=localequations(2,ioffq+ieq)
3909  rmeas=real(glder(ja),mpd) ! data
3910  rerr =real(glder(jb),mpd) ! ... and the error
3911  wght =1.0_mpd/rerr**2 ! weight from error
3912  nweig=nweig+1
3913  resid=rmeas-localcorrections(ioffq+ieq) ! subtract previous fit
3914  IF(nloopn /= 1.AND.iter /= 1.AND.lhuber /= 0) THEN
3915  IF(iter <= 3) THEN
3916  IF(abs(resid) > chuber*rerr) THEN ! down-weighting
3917  wght=wght*chuber*rerr/abs(resid)
3918  ndown=ndown+1
3919  END IF
3920  ELSE ! Cauchy
3921  wght=wght/(1.0+(resid/rerr/cauchy)**2)
3922  END IF
3923  END IF
3924 
3925  IF(lprnt.AND.iter /= 1.AND.nter /= 1) THEN
3926  chast=' '
3927  IF(abs(resid) > chuber*rerr) chast='* '
3928  IF(abs(resid) > 3.0*rerr) chast='** '
3929  IF(abs(resid) > 6.0*rerr) chast='***'
3930  IF(imeas == 0) WRITE(1,*) 'Second loop: accumulate'
3931  IF(imeas == 0) WRITE(1,103)
3932  imeas=imeas+1
3933  down=1.0/sqrt(wght)
3934  r1=resid/rerr
3935  r2=resid/down
3936  WRITE(1,104) imeas,rmeas,resid,rerr,r1,chast,r2
3937  END IF
3938 103 FORMAT(' index corrvalue residuum sigma', &
3939  ' nresid cnresid')
3940 104 FORMAT(i6,2x,2g12.4,' +-',g12.4,f7.2,1x,a3,f8.2)
3941 
3942  DO j=1,jb-ja-1 ! normal equations, local parameter loop
3943  ij=inder(ja+j) ! local parameter index J
3944  blvec(ij)=blvec(ij)+wght*rmeas*real(glder(ja+j),mpd)
3945  DO k=1,j
3946  ik=inder(ja+k) ! local parameter index K
3947  jk=ijsym(ij,ik) ! index in symmetric matrix
3948  clmat(jk)=clmat(jk) & ! force double precision
3949  +wght*real(glder(ja+j),mpd)*real(glder(ja+k),mpd)
3950  ! check for band matrix substructure
3951  IF (iter == 1) THEN
3952  id=iabs(ij-ik)+1
3953  im=min(ij,ik) ! upper/left border
3954  ibandh(id)=max(ibandh(id),im)
3955  im=min(nalc+1-ij,nalc+1-ik) ! lower/rght border (mirrored)
3956  ibandh(nalc+id)=max(ibandh(nalc+id),im)
3957  END IF
3958  END DO
3959  END DO
3960  END DO
3961  ! for non trivial fits check for bordered band matrix structure
3962  IF (iter == 1.AND.nalc > 5.AND.lfitbb > 0) THEN
3963  kx=-1
3964  kbdrx=0
3965  icmn=int(nalc,mpl)**3 ! cost (*6) should improve by at least factor 2
3966  ! upper/left border ?
3967  kbdr=0
3968  DO k=nalc,2,-1
3969  kbnd=k-2
3970  kbdr=max(kbdr,ibandh(k))
3971  icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
3972  IF (icost < icmn) THEN
3973  icmn=icost
3974  kx=k
3975  kbdrx=kbdr
3976  mside=1
3977  END IF
3978  END DO
3979  IF (kx < 0) THEN
3980  ! lower/right border instead?
3981  kbdr=0
3982  DO k=nalc,2,-1
3983  kbnd=k-2
3984  kbdr=max(kbdr,ibandh(k+nalc))
3985  icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
3986  IF (icost < icmn) THEN
3987  icmn=icost
3988  kx=k
3989  kbdrx=kbdr
3990  mside=2
3991  END IF
3992  END DO
3993  END IF
3994  IF (kx > 0) THEN
3995  mbnd=kx-2
3996  mbdr=kbdrx
3997  END IF
3998  END IF
3999 
4000  IF (mbnd >= 0) THEN
4001  ! fast solution for border banded matrix (inverse for ICALCM>0)
4002  IF (nloopn == 1) THEN
4003  nbndr(mside)=nbndr(mside)+1
4004  nbdrx=max(nbdrx,mbdr)
4005  nbndx=max(nbndx,mbnd)
4006  END IF
4007 
4008  inv=0
4009  IF (nloopn <= lfitnp.AND.iter == 1) inv=1 ! band part of inverse (for pulls)
4010  IF (icalcm == 1.OR.lprnt) inv=2 ! complete inverse
4011  IF (mside == 1) THEN
4012  CALL sqmibb(clmat,blvec,nalc,mbdr,mbnd,inv,nrank, &
4014  ELSE
4015  CALL sqmibb2(clmat,blvec,nalc,mbdr,mbnd,inv,nrank, &
4017  ENDIF
4018  ELSE
4019  ! full inversion and solution
4020  inv=2
4021  CALL sqminv(clmat,blvec,nalc,nrank,scdiag,scflag)
4022  END IF
4023  ! check for NaNs
4024  nan=0
4025  DO k=1, nalc
4026  IF ((.NOT.(blvec(k) <= 0.0_mpd)).AND. (.NOT.(blvec(k) > 0.0_mpd))) nan=nan+1
4027  END DO
4028 
4029  IF(lprnt) THEN
4030  WRITE(1,*) ' '
4031  WRITE(1,*) 'Parameter determination:',nalc,' parameters,', ' rank=',nrank
4032  WRITE(1,*) '-----------------------'
4033  IF(ndown /= 0) WRITE(1,*) ' ',ndown,' data down-weighted'
4034  WRITE(1,*) ' '
4035  END IF
4036 
4037  ! ----- third loop -------------------------------------------------
4038  ! calculate single residuals remaining after local fit and chi^2
4039 
4040  summ=0.0_mpd
4041  suwt=0.0
4042  imeas=0
4043  DO ieq=1,neq! loop over measurements
4044  ja=localequations(1,ioffq+ieq)
4045  jb=localequations(2,ioffq+ieq)
4046  ist=localequations(3,ioffq+ieq)
4047  rmeas=real(glder(ja),mpd) ! data (global contrib. subtracted)
4048  rerr =real(glder(jb),mpd) ! ... and the error
4049  wght =1.0_mpd/rerr**2 ! weight from error
4050  rmloc=0.0 ! local fit result reset
4051  DO j=1,jb-ja-1 ! local parameter loop
4052  ij=inder(ja+j)
4053  rmloc=rmloc+real(glder(ja+j),mpd)*blvec(ij) ! local fit result
4054  END DO
4055  localcorrections(ioffq+ieq)=rmloc ! save local fit result
4056  rmeas=rmeas-rmloc ! reduced to residual
4057 
4058  ! calculate pulls? (needs covariance matrix)
4059  IF(iter == 1.AND.inv > 0.AND.nloopn <= lfitnp) THEN
4060  dvar=0.0_mpd
4061  DO j=1,jb-ja-1
4062  ij=inder(ja+j)
4063  DO k=1,jb-ja-1
4064  ik=inder(ja+k)
4065  jk=ijsym(ij,ik)
4066  dvar=dvar+clmat(jk)*real(glder(ja+j),mpd)*real(glder(ja+k),mpd)
4067  END DO
4068  END DO
4069  ! some variance left to define a pull?
4070  IF (0.999999_mpd/wght > dvar) THEN
4071  pull=rmeas/sqrt(1.0_mpd/wght-dvar)
4072  IF (lhist) THEN
4073  IF (jb < ist) THEN
4074  CALL hmpent(13,real(pull,mps)) ! histogram pull
4075  CALL gmpms(5,rec,real(pull,mps))
4076  ELSE
4077  CALL hmpent(14,real(pull,mps)) ! histogram pull
4078  END IF
4079  END IF
4080  ! monitoring
4081  IF (imonit /= 0) THEN
4082  IF (jb < ist) THEN
4083  ij=inder(jb+1) ! group by first global label
4084  if (imonmd == 0) THEN
4085  irbin=min(measbins,max(1,int(pull*rerr/measres(ij)/measbinsize+0.5*real(measbins,mpd))))
4086  ELSE
4087  irbin=min(measbins,max(1,int(pull/measbinsize+0.5*real(measbins,mpd))))
4088  ENDIF
4089  irbin=irbin+measbins*(measindex(ij)-1+nummeas*iproc)
4090  meashists(irbin)=meashists(irbin)+1
4091  ENDIF
4092  ENDIF
4093  END IF
4094  END IF
4095 
4096  IF(iter == 1.AND.jb < ist.AND.lhist) &
4097  CALL gmpms(4,rec,real(rmeas/rerr,mps)) ! residual (with global deriv.)
4098 
4099  dchi2=wght*rmeas*rmeas
4100  ! DCHIT=DCHI2
4101  resid=rmeas
4102  IF(nloopn /= 1.AND.iter /= 1.AND.lhuber /= 0) THEN
4103  IF(iter <= 3) THEN
4104  IF(abs(resid) > chuber*rerr) THEN ! down-weighting
4105  wght=wght*chuber*rerr/abs(resid)
4106  dchi2=2.0*chuber*(abs(resid)/rerr-0.5*chuber)
4107  END IF
4108  ELSE
4109  wght=wght/(1.0_mpd+(resid/rerr/cauchy)**2)
4110  dchi2=log(1.0_mpd+(resid/rerr/cauchy)**2)*cauchy**2
4111  END IF
4112  END IF
4113 
4114  down=1.0/sqrt(wght)
4115 
4116  ! SUWT=SUWT+DCHI2/DCHIT
4117  suwt=suwt+rerr/down
4118  IF(lprnt) THEN
4119  chast=' '
4120  IF(abs(resid) > chuber*rerr) chast='* '
4121  IF(abs(resid) > 3.0*rerr) chast='** '
4122  IF(abs(resid) > 6.0*rerr) chast='***'
4123  IF(imeas == 0) WRITE(1,*) 'Third loop: single residuals'
4124  IF(imeas == 0) WRITE(1,105)
4125  imeas=imeas+1
4126  r1=resid/rerr
4127  r2=resid/down
4128  IF(resid < 0.0) r1=-r1
4129  IF(resid < 0.0) r2=-r2
4130  WRITE(1,106) imeas,glder(ja),rmeas,rerr,r1,chast,r2
4131  END IF
4132 105 FORMAT(' index corrvalue residuum sigma', &
4133  ' nresid cnresid')
4134 106 FORMAT(i6,2x,2g12.4,' +-',g12.4,f7.2,1x,a3,f8.2)
4135 
4136  IF(iter == nter) THEN
4137  readbufferdatad(ja)=real(rmeas,mpr8) ! store remaining residual
4138  resmax=max(resmax,abs(rmeas)/rerr)
4139  END IF
4140 
4141  IF(iter == 1.AND.lhist) THEN
4142  IF (jb < ist) THEN
4143  CALL hmpent( 3,real(rmeas/rerr,mps)) ! histogram norm residual
4144  ELSE
4145  CALL hmpent(12,real(rmeas/rerr,mps)) ! histogram norm residual
4146  END IF
4147  END IF
4148  summ=summ+dchi2 ! accumulate chi-square sum
4149  END DO
4150 
4151  ndf=neq-nrank
4152  resing=(real(nweig,mps)-real(suwt,mps))/real(nweig,mps)
4153  IF (lhist) THEN
4154  IF(iter == 1) CALL hmpent( 5,real(ndf,mps)) ! histogram Ndf
4155  IF(iter == 1) CALL hmpent(11,real(nalc,mps)) ! histogram Nlocal
4156  IF(nloopn == 2.AND.iter == nter) CALL hmpent(6,resing)
4157  END IF
4158  IF(lprnt) THEN
4159  WRITE(1,*) ' '
4160  WRITE(1,*) 'Chi^2=',summ,' at',ndf,' degrees of freedom: ', &
4161  '3-sigma limit is',chindl(3,ndf)*real(ndf,mps)
4162  WRITE(1,*) suwt,' is sum of factors, compared to',nweig, &
4163  ' Downweight fraction:',resing
4164  END IF
4165  IF(nrank /= nalc.OR.nan > 0) THEN
4166  nrej(0)=nrej(0)+1 ! count cases
4167  IF (nrec3 == huge(nrec3)) nrec3=nrc
4168  IF(lprnt) THEN
4169  WRITE(1,*) ' rank deficit/NaN ', nalc, nrank, nan
4170  WRITE(1,*) ' ---> rejected!'
4171  END IF
4172  GO TO 90
4173  END IF
4174  IF(ndf <= 0) THEN
4175  nrej(1)=nrej(1)+1 ! count cases
4176  IF(lprnt) THEN
4177  WRITE(1,*) ' ---> rejected!'
4178  END IF
4179  GO TO 90
4180  END IF
4181 
4182  chndf=real(summ/real(ndf,mpd),mps)
4183 
4184  IF(iter == 1.AND.lhist) CALL hmpent(4,chndf) ! histogram chi^2/Ndf
4185  END DO ! outlier iteration loop
4186 
4187  ndfs=ndfs+ndf ! (local) sum of Ndf
4188  sndf=sndf+real(ndf,mpd)*dw1 ! (local) weighted sum of Ndf
4189 
4190  ! ----- reject eventually ------------------------------------------
4191 
4192  IF(newite.AND.iterat == 2) THEN ! find record with largest Chi^2/Ndf
4193  IF(nrecp2 < 0.AND.chndf > writebufferdata(2,iproc+1)) THEN
4194  writebufferdata(2,iproc+1)=chndf
4195  writebufferinfo(7,iproc+1)=nrc
4196  END IF
4197  END IF
4198 
4199  chichi=chindl(3,ndf)*real(ndf,mps)
4200  ! GF IF(SUMM.GT.50.0*CHICHI) THEN ! huge
4201  ! CHK CHICUT<0: NO cut (1st iteration)
4202  IF(chicut >= 0.0) THEN
4203  IF(summ > chhuge*chichi) THEN ! huge
4204  nrej(2)=nrej(2)+1 ! count cases with huge chi^2
4205  IF(lprnt) THEN
4206  WRITE(1,*) ' ---> rejected!'
4207  END IF
4208  GO TO 90
4209  END IF
4210 
4211  IF(chicut > 0.0) THEN
4212  chlimt=chicut*chichi
4213  ! WRITE(*,*) 'chi^2 ',SUMM,CHLIMT,CHICUT,CHINDL(3,NDF),NDF
4214  IF(summ > chlimt) THEN
4215  IF(lprnt) THEN
4216  WRITE(1,*) ' ---> rejected!'
4217  END IF
4218  ! add to FVALUE
4219  dchi2=chlimt ! total contribution limit
4220  dchi2s=dchi2s+dchi2*dw1 ! add total contribution
4221  nrej(3)=nrej(3)+1 ! count cases with large chi^2
4222  GO TO 90
4223  END IF
4224  END IF
4225  END IF
4226 
4227  IF(lhuber > 1.AND.dwcut /= 0.0.AND.resing > dwcut) THEN
4228  ! add to FVALUE
4229  dchi2=summ ! total contribution
4230  dchi2s=dchi2s+dchi2*dw1 ! add total contribution
4231  nrej(3)=nrej(3)+1 ! count cases with large chi^2
4232  ! WRITE(*,*) 'Downweight fraction cut ',RESING,DWCUT,SUMM
4233  IF(lprnt) THEN
4234  WRITE(1,*) ' ---> rejected!'
4235  END IF
4236  GO TO 90
4237  END IF
4238 
4239  IF(newite.AND.iterat == 2) THEN ! find record with largest residual
4240  IF(nrecpr < 0.AND.resmax > writebufferdata(1,iproc+1)) THEN
4241  writebufferdata(1,iproc+1)=real(resmax,mps)
4242  writebufferinfo(6,iproc+1)=nrc
4243  END IF
4244  END IF
4245  ! 'track quality' per binary file: accepted records
4246  naccf(kfl)=naccf(kfl)+1
4247  ndff(kfl) =ndff(kfl) +ndf
4248  chi2f(kfl)=chi2f(kfl)+chndf
4249 
4250  ! ----- fourth loop ------------------------------------------------
4251  ! update of global matrix and vector according to the "Millepede"
4252  ! principle, from the global/local information
4253 
4254  DO ieq=1,neq! loop over measurements
4255  ja=localequations(1,ioffq+ieq)
4256  jb=localequations(2,ioffq+ieq)
4257  ist=localequations(3,ioffq+ieq)
4258  rmeas=real(glder(ja),mpd) ! data residual
4259  rerr =real(glder(jb),mpd) ! ... and the error
4260  wght =1.0_mpd/rerr**2 ! weight from measurement error
4261  dchi2=wght*rmeas*rmeas ! least-square contribution
4262 
4263  IF(nloopn /= 1.AND.lhuber /= 0) THEN ! check residual
4264  resid=abs(rmeas)
4265  IF(resid > chuber*rerr) THEN
4266  wght=wght*chuber*rerr/resid ! down-weighting
4267  dchi2=2.0*chuber*(resid/rerr-0.5*chuber) ! modified contribution
4268  END IF
4269  END IF
4270  dchi2s=dchi2s+dchi2*dw1 ! add to total objective function
4271 
4272  ! global-global matrix contribution: add directly to gg-matrix
4273 
4274  DO j=1,ist-jb
4275  ivgbj=globalparlabelindex(2,inder(jb+j)) ! variable-parameter index
4276  IF(ivgbj > 0) THEN
4277  globalvector(ioffb+ivgbj)=globalvector(ioffb+ivgbj) &
4278  +dw1*wght*rmeas*real(glder(jb+j),mpd) ! vector !!! reverse
4279  globalcounter(ioffb+ivgbj)=globalcounter(ioffb+ivgbj)+1
4280  IF(icalcm == 1) THEN
4281  ije=backindexusage(ioffe+ivgbj) ! get index of index, non-zero
4282  DO k=1,j
4283  ivgbk=globalparlabelindex(2,inder(jb+k))
4284  IF(ivgbk > 0) THEN
4285  ike=backindexusage(ioffe+ivgbk) ! get index of index, non-zero
4286  ia=max(ije,ike) ! larger
4287  ib=min(ije,ike) ! smaller
4288  ij=ib+(ia*ia-ia)/2
4289  writebufferupdates(ioffd+ij)=writebufferupdates(ioffd+ij) &
4290  -dw1*wght*real(glder(jb+j),mpd)*real(glder(jb+k),mpd)
4291  END IF
4292  END DO
4293  END IF
4294  END IF
4295  END DO
4296 
4297  ! normal equations - rectangular matrix for global/local pars
4298  ! global-local matrix contribution: accumulate rectangular matrix
4299  IF (icalcm /= 1) cycle
4300  DO j=1,ist-jb
4301  ivgbj=globalparlabelindex(2,inder(jb+j)) ! variable-parameter index
4302  IF(ivgbj > 0) THEN
4303  ije=backindexusage(ioffe+ivgbj) ! get index of index, non-zero
4304  DO k=1,jb-ja-1
4305  ik=inder(ja+k) ! local index
4306  jk=ik+(ije-1)*nalc ! matrix index
4307  localglobalmatrix(jk)=localglobalmatrix(jk)+dw2*wght*real(glder(jb+j),mpd)*real(glder(ja+k),mpd)
4308  localglobalmap(jk)=localglobalmap(jk)+1
4309  END DO
4310  END IF
4311  END DO
4312  END DO
4313 
4314  ! ----- final matrix update ----------------------------------------
4315  ! update global matrices and vectors
4316  IF(icalcm /= 1) GO TO 90 ! matrix update
4317  ! (inverse local matrix) * (rectang. matrix) -> CORM
4318  ! T
4319  ! resulting symmetrix matrix = G * Gamma^{-1} * G
4320 
4321  ! check sparsity of localGlobalMatrix (with par. groups)
4322  isize=nalc+nalg+1 ! row/clolumn offsets
4323  ! check rows
4324  k=0 ! offset
4325  DO i=1, nalg
4326  localglobalstructure(i)=isize
4327  DO j=1, nalc
4328  IF (localglobalmap(k+j) > 0) THEN
4329  localglobalstructure(isize+1)=j ! column
4330  localglobalstructure(isize+2)=k+j ! index
4331  isize=isize+2
4332  ENDIF
4333  END DO
4334  k=k+nalc
4335  END DO
4336  ! <50% non-zero elements?
4337  IF (isize-localglobalstructure(1) < nalc*nalg) THEN
4338  ! check columns (too)
4339  DO j=1, nalc
4340  localglobalstructure(nalg+j)=isize
4341  k=0 ! offset
4342  DO i=1, nalg
4343  IF (localglobalmap(k+j) > 0) THEN
4344  localglobalstructure(isize+1)=i ! row
4345  localglobalstructure(isize+2)=k+j ! index
4346  isize=isize+2
4347  ENDIF
4348  k=k+nalc
4349  END DO
4350  END DO
4351  localglobalstructure(nalg+nalc+1)=isize
4353  ELSE
4354  CALL dbavat(clmat,localglobalmatrix,writebufferupdates(ioffd+1),nalc,-nalg)
4355  END IF
4356  ! (rectang. matrix) * (local param vector) -> CORV
4357  ! resulting vector = G * q (q = local parameter)
4358  ! CALL DBGAX(DQ(IGLMA/2+1),BLVEC,DQ(ICORV/2+1),NALG,NALC) ! not done
4359  ! the vector update is not done, because after local fit it is zero!
4360 
4361  ! update cache status
4362  writebufferinfo(1,iproc+1)=writebufferinfo(1,iproc+1)+1
4363  writebufferinfo(2,iproc+1)=writebufferinfo(2,iproc+1)+ngg
4364  writebufferinfo(3,iproc+1)=writebufferinfo(3,iproc+1)+ngrp+2
4365  ! check free space
4366  nfred=writebufferheader(-1)-writebufferinfo(2,iproc+1)-writebufferheader(-2)
4367  nfrei=writebufferheader(1)-writebufferinfo(3,iproc+1)-writebufferheader(2)
4368  IF (nfred < 0.OR.nfrei < 0) THEN ! need to flush
4369  nb=writebufferinfo(1,iproc+1)
4370  joffd=writebufferheader(-1)*iproc ! offset data
4371  joffi=writebufferheader(1)*iproc+2 ! offset indices
4372  used=real(writebufferinfo(2,iproc+1),mps)/real(writebufferheader(-1),mps)
4373  writebufferinfo(4,iproc+1)=writebufferinfo(4,iproc+1) +nint(1000.0*used,mpi)
4374  used=real(writebufferinfo(3,iproc+1),mps)/real(writebufferheader(1),mps)
4375  writebufferinfo(5,iproc+1)=writebufferinfo(5,iproc+1) +nint(1000.0*used,mpi)
4376  !$OMP CRITICAL
4379 
4380  DO ib=1,nb
4381  il=1 ! row in update matrix
4382  DO in=1,writebufferindices(joffi)
4383  i=writebufferindices(joffi+in)
4384  j=writebufferindices(joffi+1) ! 1. group
4385  iprc=ijprec(i,j) ! group pair precision
4386  jl=1 ! col in update matrix
4387  ! start (rows) for continous groups
4388  j1=j
4389  jl1=jl
4390  ! other groups for row
4391  DO jn=2,in
4393  jnx=writebufferindices(joffi+jn) ! next group
4394  iprcnx=ijprec(i,jnx) ! group pair precision
4395  ! end of continous groups?
4396  IF (.NOT.((jnx == j+1).AND.(iprc == iprcnx))) THEN
4397  CALL mgupdt(i,j1,j,il,jl1,writebufferupdates(joffd+1)) ! matrix update
4398  !print *, ' update ', ib,i,j1,j,il,jl1,0,iprc,jnx,iprcnx
4399  ! restart continous groups
4400  j1=jnx ! new 1. column
4401  jl1=jl
4402  iprc=iprcnx
4403  END IF
4404  j=jnx ! last group
4405  END DO
4406  CALL mgupdt(i,j1,j,il,jl1,writebufferupdates(joffd+1)) ! final matrix update
4407  !print *, '.update ', ib, i,j1,j,il,jl1,1,iprc
4409  END DO
4410  joffd=joffd+(il*il-il)/2
4411  joffi=joffi+writebufferindices(joffi)+2
4412  END DO
4413  !$OMP END CRITICAL
4414  ! reset counter, pointers
4415  DO k=1,3
4416  writebufferinfo(k,iproc+1)=0
4417  END DO
4418  END IF
4419 
4420 90 IF(lprnt) THEN
4421  WRITE(1,*) ' '
4422  WRITE(1,*) '------------------ End of printout for record',nrc
4423  WRITE(1,*) ' '
4424  END IF
4425 
4426  DO i=1,nalg ! reset global index array
4427  iext=globalindexusage(ioffc+i)
4428  backindexusage(ioffe+iext)=0
4429  END DO
4430 
4431  END DO
4432  !$OMP END PARALLEL DO
4433 
4434  IF (icalcm == 1) THEN
4435  ! flush remaining matrices
4436  DO k=1,mthrd ! update statistics
4438  used=real(writebufferinfo(2,k),mps)/real(writebufferheader(-1),mps)
4439  writebufferinfo(4,k)=writebufferinfo(4,k)+nint(1000.0*used,mpi)
4442  writebufferinfo(4,k)=0
4444  used=real(writebufferinfo(3,k),mps)/real(writebufferheader(1),mps)
4445  writebufferinfo(5,k)=writebufferinfo(5,k)+nint(1000.0*used,mpi)
4448  writebufferinfo(5,k)=0
4449  END DO
4450 
4451  !$OMP PARALLEL &
4452  !$OMP DEFAULT(PRIVATE) &
4453  !$OMP SHARED(writeBufferHeader,writeBufferInfo,writeBufferIndices,writeBufferUpdates,MTHRD) &
4454  !$OMP SHARED(globalAllParToGroup,globalAllIndexGroups,nspc)
4455  iproc=0
4456  !$ IPROC=OMP_GET_THREAD_NUM() ! thread number
4457  DO jproc=0,mthrd-1
4458  nb=writebufferinfo(1,jproc+1)
4459  ! print *, ' flush end ', JPROC, NRC, NB
4460  joffd=writebufferheader(-1)*jproc ! offset data
4461  joffi=writebufferheader(1)*jproc+2 ! offset indices
4462  DO ib=1,nb
4463  ! print *, ' buf end ', JPROC,IB,writeBufferIndices(JOFFI-1),writeBufferIndices(JOFFI)
4464  il=1 ! row in update matrix
4465  DO in=1,writebufferindices(joffi)
4466  i=writebufferindices(joffi+in)
4467  !$ IF (MOD(I,MTHRD).EQ.IPROC) THEN
4468  j=writebufferindices(joffi+1) ! 1. group
4469  iprc=ijprec(i,j) ! group pair precision
4470  jl=1 ! col in update matrix
4471  ! start (rows) for continous groups
4472  j1=j
4473  jl1=jl
4474  ! other groups for row
4475  DO jn=2,in
4477  jnx=writebufferindices(joffi+jn) ! next group
4478  iprcnx=ijprec(i,jnx) ! group pair precision
4479  ! end of continous groups?
4480  IF (.NOT.((jnx == j+1).AND.(iprc == iprcnx))) THEN
4481  CALL mgupdt(i,j1,j,il,jl1,writebufferupdates(joffd+1)) ! matrix update
4482  !print *, ' update ', ib,i,j1,j,il,jl1,0,iprc,jnx,iprcnx
4483  ! restart continous groups
4484  j1=jnx ! new 1. column
4485  jl1=jl
4486  iprc=iprcnx
4487  END IF
4488  j=jnx ! last group
4489  END DO
4490  CALL mgupdt(i,j1,j,il,jl1,writebufferupdates(joffd+1)) ! final matrix update
4491  !print *, '.update ', ib, i,j1,j,il,jl1,1,iprc
4492  !$ END IF
4494  END DO
4495  joffd=joffd+(il*il-il)/2
4496  joffi=joffi+writebufferindices(joffi)+2
4497  END DO
4498  END DO
4499  !$OMP END PARALLEL
4500  END IF
4501 
4502  IF(newite.AND.iterat == 2) THEN ! get worst records (for printrecord -1 -1)
4503  IF (nrecpr < 0) THEN
4504  DO k=1,mthrd
4505  IF (writebufferdata(1,k) > value1) THEN
4506  value1=writebufferdata(1,k)
4507  nrec1 =writebufferinfo(6,k)
4508  END IF
4509  END DO
4510  END IF
4511  IF (nrecp2 < 0) THEN
4512  DO k=1,mthrd
4513  IF (writebufferdata(2,k) > value2) THEN
4514  value2=writebufferdata(2,k)
4515  nrec2 =writebufferinfo(7,k)
4516  END IF
4517  END DO
4518  END IF
4519  END IF
4520 
4521 END SUBROUTINE loopbf
4522 
4523 
4524 
4525 
4526 !***********************************************************************
4527 
4540 SUBROUTINE prtglo
4541  USE mpmod
4542 
4543  IMPLICIT NONE
4544  REAL(mps):: dpa
4545  REAL(mps):: err
4546  REAL(mps):: gcor
4547  INTEGER(mpi) :: i
4548  INTEGER(mpi) :: icount
4549  INTEGER(mpi) :: ie
4550  INTEGER(mpi) :: iev
4551  INTEGER(mpi) :: ij
4552  INTEGER(mpi) :: imin
4553  INTEGER(mpi) :: iprlim
4554  INTEGER(mpi) :: isub
4555  INTEGER(mpi) :: itgbi
4556  INTEGER(mpi) :: itgbl
4557  INTEGER(mpi) :: ivgbi
4558  INTEGER(mpi) :: j
4559  INTEGER(mpi) :: label
4560  INTEGER(mpi) :: lup
4561  REAL(mps):: par
4562  LOGICAL :: lowstat
4563 
4564  REAL(mpd):: diag
4565  REAL(mpd)::gmati
4566  REAL(mpd)::gcor2
4567  INTEGER(mpi) :: labele(3)
4568  REAL(mps):: compnt(3)
4569  SAVE
4570  ! ...
4571 
4572  lup=09
4573  CALL mvopen(lup,'millepede.res')
4574 
4575  WRITE(*,*) ' '
4576  WRITE(*,*) ' Result of fit for global parameters'
4577  WRITE(*,*) ' ==================================='
4578  WRITE(*,*) ' '
4579 
4580  WRITE(*,101)
4581 
4582  WRITE(lup,*) 'Parameter ! first 3 elements per line are', &
4583  ' significant (if used as input)'
4584 
4585 
4586  iprlim=10
4587  DO itgbi=1,ntgb ! all parameter variables
4588  itgbl=globalparlabelindex(1,itgbi)
4589  ivgbi=globalparlabelindex(2,itgbi)
4590  par=real(globalparameter(itgbi),mps) ! initial value
4591  icount=0 ! counts
4592  lowstat = .false.
4593  IF(ivgbi > 0) THEN
4594  icount=globalcounter(ivgbi) ! used in last iteration
4595  lowstat = (icount < mreqena) ! to few accepted entries
4596  dpa=real(globalparameter(itgbi)-globalparstart(itgbi),mps) ! difference
4597  IF(ALLOCATED(workspacediag)) THEN ! provide parameter errors?
4598  gmati=globalmatd(globalrowoffsets(ivgbi)+ivgbi)
4599  err=sqrt(abs(real(gmati,mps)))
4600  IF(gmati < 0.0_mpd) err=-err
4601  diag=workspacediag(ivgbi)
4602  gcor=-1.0
4603  IF(gmati*diag > 0.0_mpd) THEN ! global correlation
4604  gcor2=1.0_mpd-1.0_mpd/(gmati*diag)
4605  IF(gcor2 >= 0.0_mpd.AND.gcor2 <= 1.0_mpd) gcor=real(sqrt(gcor2),mps)
4606  END IF
4607  END IF
4608  END IF
4609  IF(ipcntr > 1) icount=globalparcounts(itgbi) ! from binary files
4610  IF(lowstat) icount=-(icount+1) ! flag 'lowstat' with icount < 0
4611  IF(itgbi <= iprlim) THEN
4612  IF(ivgbi <= 0) THEN
4613  WRITE(* ,102) itgbl,par,real(globalparpresigma(itgbi),mps)
4614  ELSE
4615  IF(ALLOCATED(workspacediag)) THEN ! provide parameter errors?
4616  IF (igcorr == 0) THEN
4617  WRITE(*,102) itgbl,par,real(globalparpresigma(itgbi),mps),dpa,err
4618  ELSE
4619  WRITE(*,102) itgbl,par,real(globalparpresigma(itgbi),mps),dpa,err,gcor
4620  END IF
4621  ELSE
4622  WRITE(*,102) itgbl,par,real(globalparpresigma(itgbi),mps),dpa
4623  END IF
4624  END IF
4625  ELSE IF(itgbi == iprlim+1) THEN
4626  WRITE(* ,*) '... (further printout suppressed, but see log file)'
4627  END IF
4628 
4629  ! file output
4630  IF(ivgbi <= 0) THEN
4631  IF (ipcntr /= 0) THEN
4632  WRITE(lup,110) itgbl,par,real(globalparpresigma(itgbi),mps),icount
4633  ELSE
4634  WRITE(lup,102) itgbl,par,real(globalparpresigma(itgbi),mps)
4635  END IF
4636  ELSE
4637  IF(ALLOCATED(workspacediag)) THEN ! provide parameter errors?
4638  IF (ipcntr /= 0) THEN
4639  WRITE(lup,112) itgbl,par,real(globalparpresigma(itgbi),mps),dpa,err,icount
4640  ELSE IF (igcorr /= 0) THEN
4641  WRITE(lup,102) itgbl,par,real(globalparpresigma(itgbi),mps),dpa,err,gcor
4642  ELSE
4643  WRITE(lup,102) itgbl,par,real(globalparpresigma(itgbi),mps),dpa,err
4644  END IF
4645  ELSE
4646  IF (ipcntr /= 0) THEN
4647  WRITE(lup,111) itgbl,par,real(globalparpresigma(itgbi),mps),dpa,icount
4648  ELSE
4649  WRITE(lup,102) itgbl,par,real(globalparpresigma(itgbi),mps),dpa
4650  END IF
4651  END IF
4652  END IF
4653  END DO
4654  rewind lup
4655  CLOSE(unit=lup)
4656 
4657  IF(metsol == 2) THEN ! diagonalisation: write eigenvectors
4658  CALL mvopen(lup,'millepede.eve')
4659  imin=1
4660  DO i=nagb,1,-1
4661  IF(workspaceeigenvalues(i) > 0.0_mpd) THEN
4662  imin=i ! index of smallest pos. eigenvalue
4663  EXIT
4664  ENDIF
4665  END DO
4666  iev=0
4667 
4668  DO isub=0,min(15,imin-1)
4669  IF(isub < 10) THEN
4670  i=imin-isub
4671  ELSE
4672  i=isub-9
4673  END IF
4674 
4675  ! DO I=IMIN,MAX(1,IMIN-9),-1 ! backward loop, up to 10 vectors
4676  WRITE(*,*) 'Eigenvector ',i,' with eigenvalue',workspaceeigenvalues(i)
4677  WRITE(lup,*) 'Eigenvector ',i,' with eigenvalue',workspaceeigenvalues(i)
4678  DO j=1,nagb
4679  ij=j+(i-1)*nagb ! index with eigenvector array
4680  IF(j <= nvgb) THEN
4681  itgbi=globalparvartototal(j)
4682  label=globalparlabelindex(1,itgbi)
4683  ELSE
4684  label=nvgb-j ! label negative for constraints
4685  END IF
4686  iev=iev+1
4687  labele(iev)=label
4688  compnt(iev)=real(workspaceeigenvectors(ij),mps) ! component
4689  IF(iev == 3) THEN
4690  WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
4691  iev=0
4692  END IF
4693  END DO
4694  IF(iev /= 0) WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
4695  iev=0
4696  WRITE(lup,*) ' '
4697  END DO
4698 
4699  END IF
4700 
4701 101 FORMAT(1x,' label parameter presigma differ', &
4702  ' error'/ 1x,'-----------',4x,4('-------------'))
4703 102 FORMAT(i10,2x,4g14.5,f8.3)
4704 103 FORMAT(3(i11,f11.7,2x))
4705 110 FORMAT(i10,2x,2g14.5,28x,i12)
4706 111 FORMAT(i10,2x,3g14.5,14x,i12)
4707 112 FORMAT(i10,2x,4g14.5,i12)
4708 END SUBROUTINE prtglo ! print final log file
4709 
4710 !***********************************************************************
4711 
4721 SUBROUTINE prtstat
4722  USE mpmod
4723 
4724  IMPLICIT NONE
4725  REAL(mps):: par
4726  REAL(mps):: presig
4727  INTEGER(mpi) :: icount
4728  INTEGER(mpi) :: itgbi
4729  INTEGER(mpi) :: itgbl
4730  INTEGER(mpi) :: itpgrp
4731  INTEGER(mpi) :: ivgbi
4732  INTEGER(mpi) :: lup
4733  INTEGER(mpi) :: ncon
4734  INTEGER(mpi) :: k
4735  CHARACTER :: c1
4736 
4737  SAVE
4738  ! ...
4739 
4740  lup=09
4741  CALL mvopen(lup,'millepede.res')
4742  WRITE(lup,*) '*** Results of checking input only, no solution performed ***'
4743  WRITE(lup,*) '! fixed-1: by pre-sigma, -2: by entries cut, -3: by iterated entries cut'
4744  WRITE(lup,*) '! Label Value Pre-sigma Entries Constraints Status '
4745  !iprlim=10
4746  DO itgbi=1,ntgb ! all parameter variables
4747  itgbl=globalparlabelindex(1,itgbi)
4748  ivgbi=globalparlabelindex(2,itgbi)
4749  c1=' '
4750  IF (globalparlabelindex(3,itgbi) == itgbl) c1='>'
4751  par=real(globalparameter(itgbi),mps) ! initial value
4752  presig=real(globalparpresigma(itgbi),mps) ! initial presigma
4753  icount=globalparcounts(itgbi) ! from binary files
4754  ncon=globalparcons(itgbi) ! number of active constraints
4755 
4756  IF (ivgbi <= 0) THEN
4757  WRITE(lup,110) c1,itgbl,par,presig,icount,ncon,ivgbi
4758  ELSE
4759  WRITE(lup,111) c1,itgbl,par,presig,icount,ncon
4760  END IF
4761  END DO
4762  ! appearance statistics
4763  IF (icheck > 1) THEN
4764  WRITE(lup,*) '! '
4765  WRITE(lup,*) '! Appearance statistics '
4766  WRITE(lup,*) '! Label First file and record Last file and record #files #paired-par'
4767  DO itgbi=1,ntgb
4768  itpgrp=globalparlabelindex(4,itgbi)
4769  WRITE(lup,112) globalparlabelindex(1,itgbi), (appearancecounter(itgbi*5+k), k=-4,0), paircounter(itpgrp)
4770  END DO
4771  END IF
4772  rewind lup
4773  CLOSE(unit=lup)
4774 
4775 110 FORMAT(' !',a1,i10,2x,2g14.5,2i12,' fixed',i2)
4776 111 FORMAT(' !',a1,i10,2x,2g14.5,2i12,' variable')
4777 112 FORMAT(' !.',i10,6i11)
4778 END SUBROUTINE prtstat ! print input statistics
4779 
4780 
4790 
4791 SUBROUTINE avprds(n,x,b,is)
4792  USE mpmod
4793 
4794  IMPLICIT NONE
4795  INTEGER(mpi) :: i
4796  INTEGER(mpi) :: ia
4797  INTEGER(mpi) :: ib
4798  INTEGER(mpi) :: in
4799  INTEGER(mpi) :: ipg
4800  INTEGER(mpi) :: iproc
4801  INTEGER(mpi) :: ir
4802  INTEGER(mpi) :: irgn
4803  INTEGER(mpi) :: j
4804  INTEGER(mpi) :: ja
4805  INTEGER(mpi) :: jb
4806  INTEGER(mpi) :: jn
4807  INTEGER(mpi) :: jrgn
4808  INTEGER(mpi) :: lj
4809 
4810  INTEGER(mpi), INTENT(IN) :: n
4811  REAL(mpd), INTENT(IN) :: x(n)
4812  REAL(mpd), INTENT(OUT) :: b(n)
4813  INTEGER(mpi), INTENT(IN) :: is(*)
4814  INTEGER(mpl) :: k
4815  INTEGER(mpl) :: kk
4816  INTEGER(mpl) :: ku
4817  INTEGER(mpl) :: ll
4818  INTEGER(mpl) :: indij
4819  INTEGER(mpl) :: indid
4820  INTEGER(mpl) :: ij
4821  INTEGER(mpi) :: ichunk
4822  !$ INTEGER(mpi) OMP_GET_THREAD_NUM
4823  SAVE
4824  ! ...
4825  !$ DO i=1,n
4826  !$ b(i)=0.0_mpd ! reset 'global' B()
4827  !$ END DO
4828  ichunk=min((n+mthrd-1)/mthrd/8+1,1024)
4829  IF(matsto /= 2) THEN
4830  ! full symmetric matrix
4831  ! parallelize row loop
4832  ! private copy of B(N) for each thread, combined at end, init with 0.
4833  ! slot of 1024 'I' for next idle thread
4834  !$OMP PARALLEL DO &
4835  !$OMP PRIVATE(J,IJ) &
4836  !$OMP REDUCTION(+:B) &
4837  !$OMP SCHEDULE(DYNAMIC,ichunk)
4838  DO i=1,n
4839  ij=globalrowoffsets(i)
4840  b(i)=globalmatd(ij+i)*x(i)
4841  DO j=1,i-1
4842  b(j)=b(j)+globalmatd(ij+j)*x(i)
4843  b(i)=b(i)+globalmatd(ij+j)*x(j)
4844  END DO
4845  END DO
4846  !$OMP END PARALLEL DO
4847  ELSE
4848  ! sparse, compressed matrix
4849  IF(sparsematrixoffsets(2,1) /= n) THEN
4850  CALL peend(24,'Aborted, vector/matrix size mismatch')
4851  stop 'AVPRDS: mismatched vector and matrix'
4852  END IF
4853  ! parallelize row (group) loop
4854  ! slot of 1024 'I' for next idle thread
4855  !$OMP PARALLEL DO &
4856  !$OMP PRIVATE(I,IR,K,KK,LL,KU,INDID,INDIJ,J,JN,LJ) &
4857  !$OMP PRIVATE(IA,IB,IN,JA,JB,IRGN,JRGN) &
4858  !$OMP REDUCTION(+:B) &
4859  !$OMP SCHEDULE(DYNAMIC,ichunk)
4860  DO ipg=1,napgrp
4861  iproc=0
4862  !$ IPROC=OMP_GET_THREAD_NUM() ! thread number
4863  ! row group
4864  ia=globalallindexgroups(ipg) ! first (global) row
4865  ib=globalallindexgroups(ipg+1)-1 ! last (global) row
4866  in=ib-ia+1 ! number of rows
4867  ! check x region
4868  irgn=1 ! non-zero region in x
4869  DO WHILE (ia > is(2*irgn+1).AND.irgn < is(1))
4870  irgn=irgn+1
4871  END DO
4872  !
4873  ! diagonal elements
4874  IF (ia <= is(2*irgn+1).AND.ib >= is(2*irgn)) THEN
4875  b(ia:ib)=globalmatd(ia:ib)*x(ia:ib)
4876  ELSE
4877  b(ia:ib)=0.0_mpd
4878  END IF
4879  ! off-diagonals double precision
4880  ir=ipg
4881  kk=sparsematrixoffsets(1,ir) ! offset in 'd' (column lists)
4882  ll=sparsematrixoffsets(2,ir) ! offset in 'j' (matrix)
4883  ku=sparsematrixoffsets(1,ir+1)-kk
4884  indid=kk
4885  indij=ll
4886  IF (sparsematrixcolumns(indid+1) /= 0) THEN ! no compression
4887  DO i=ia,ib
4888  IF (i <= is(2*irgn+1).AND.i >= is(2*irgn)) THEN
4889  DO k=1,ku
4890  j=sparsematrixcolumns(indid+k)
4891  b(j)=b(j)+globalmatd(indij+k)*x(i)
4892  END DO
4893  END IF
4894  jrgn=1 ! non-zero region in x
4895  DO k=1,ku
4896  j=sparsematrixcolumns(indid+k)
4897  DO WHILE (j > is(2*jrgn+1).AND.jrgn < is(1))
4898  jrgn=jrgn+1
4899  END DO
4900  IF (j <= is(2*jrgn+1).AND.j >= is(2*jrgn)) THEN
4901  b(i)=b(i)+globalmatd(indij+k)*x(j)
4902  END IF
4903  END DO
4904  indij=indij+ku
4905  END DO
4906  ELSE
4907  jrgn=1 ! non-zero region in x
4908  ! regions of continous column groups
4909  DO k=2,ku-2,2
4910  j=sparsematrixcolumns(indid+k) ! first group
4911  ja=globalallindexgroups(j) ! first (global) column
4912  lj=sparsematrixcolumns(indid+k-1) ! region offset
4913  jn=sparsematrixcolumns(indid+k+1)-lj ! number of columns
4914  jb=ja+jn-1 ! last (global) column
4915  ! check x region
4916  DO WHILE (ja > is(2*jrgn+1).AND.jrgn < is(1))
4917  jrgn=jrgn+1
4918  END DO
4919  IF (ja <= is(2*jrgn+1).AND.jb >= is(2*jrgn)) THEN
4920  lj=1 ! index (in group region)
4921  DO i=ia,ib
4922  b(i)=b(i)+dot_product(globalmatd(indij+lj:indij+lj+jn-1),x(ja:jb))
4923  lj=lj+jn
4924  END DO
4925  END IF
4926  IF (mextnd == 0.AND.ia <= is(2*irgn+1).AND.ib >= is(2*irgn)) THEN
4927  lj=1
4928  DO j=ja,jb
4929  b(j)=b(j)+dot_product(globalmatd(indij+lj:indij+jn*in:jn),x(ia:ib))
4930  lj=lj+1
4931  END DO
4932  END IF
4933  indij=indij+in*jn
4934  END DO
4935  END IF
4936  ! mixed precision
4937  IF (nspc > 1) THEN
4938  ir=ipg+napgrp+1 ! off-diagonals single precision
4939  kk=sparsematrixoffsets(1,ir) ! offset in 'd' (column lists)
4940  ll=sparsematrixoffsets(2,ir) ! offset in 'j' (matrix)
4941  ku=sparsematrixoffsets(1,ir+1)-kk
4942  indid=kk
4943  indij=ll
4944  IF (sparsematrixcolumns(indid+1) /= 0) THEN ! no compression
4945  DO i=ia,ib
4946  DO k=1,ku
4947  j=sparsematrixcolumns(indid+k)
4948  b(j)=b(j)+real(globalmatf(indij+k),mpd)*x(i)
4949  b(i)=b(i)+real(globalmatf(indij+k),mpd)*x(j)
4950  END DO
4951  indij=indij+ku
4952  END DO
4953  ELSE
4954  jrgn=1 ! non-zero region in x
4955  ! regions of continous column groups
4956  DO k=2,ku-2,2
4957  j=sparsematrixcolumns(indid+k) ! first group
4958  ja=globalallindexgroups(j) ! first (global) column
4959  lj=sparsematrixcolumns(indid+k-1) ! region offset
4960  jn=sparsematrixcolumns(indid+k+1)-lj ! number of columns
4961  jb=ja+jn-1 ! last (global) column
4962  ! check x region
4963  DO WHILE (ja > is(2*jrgn+1).AND.jrgn < is(1))
4964  jrgn=jrgn+1
4965  END DO
4966  IF (ja <= is(2*jrgn+1).AND.jb >= is(2*jrgn)) THEN
4967  lj=1 ! index (in group region)
4968  DO i=ia,ib
4969  b(i)=b(i)+dot_product(real(globalmatf(indij+lj:indij+lj+jn-1),mpd),x(ja:jb))
4970  lj=lj+jn
4971  END DO
4972  END IF
4973  IF (mextnd == 0.AND.ia <= is(2*irgn+1).AND.ib >= is(2*irgn)) THEN
4974  lj=1
4975  DO j=ja,jb
4976  b(j)=b(j)+dot_product(real(globalmatf(indij+lj:indij+jn*in:jn),mpd),x(ia:ib))
4977  lj=lj+1
4978  END DO
4979  END IF
4980  indij=indij+in*jn
4981  END DO
4982  END IF
4983  END IF
4984  END DO
4985  ENDIF
4986 
4987 END SUBROUTINE avprds
4988 
5000 
5001 SUBROUTINE avprd0(n,l,x,b)
5002  USE mpmod
5003 
5004  IMPLICIT NONE
5005  INTEGER(mpi) :: i
5006  INTEGER(mpi) :: ia
5007  INTEGER(mpi) :: ib
5008  INTEGER(mpi) :: in
5009  INTEGER(mpi) :: ipg
5010  INTEGER(mpi) :: iproc
5011  INTEGER(mpi) :: ir
5012  INTEGER(mpi) :: j
5013  INTEGER(mpi) :: ja
5014  INTEGER(mpi) :: jb
5015  INTEGER(mpi) :: jn
5016  INTEGER(mpi) :: lj
5017 
5018  INTEGER(mpi), INTENT(IN) :: n
5019  INTEGER(mpl), INTENT(IN) :: l
5020  REAL(mpd), INTENT(IN) :: x(n)
5021  REAL(mpd), INTENT(OUT) :: b(n)
5022  INTEGER(mpl) :: k
5023  INTEGER(mpl) :: kk
5024  INTEGER(mpl) :: ku
5025  INTEGER(mpl) :: ll
5026  INTEGER(mpl) :: indij
5027  INTEGER(mpl) :: indid
5028  INTEGER(mpl) :: ij
5029  INTEGER(mpi) :: ichunk
5030  !$ INTEGER(mpi) OMP_GET_THREAD_NUM
5031  SAVE
5032  ! ...
5033  !$ DO i=1,n
5034  !$ b(i)=0.0_mpd ! reset 'global' B()
5035  !$ END DO
5036  ichunk=min((n+mthrd-1)/mthrd/8+1,1024)
5037  IF(matsto /= 2) THEN
5038  ! full or unpacked (block diagonal) symmetric matrix
5039  ! parallelize row loop
5040  ! private copy of B(N) for each thread, combined at end, init with 0.
5041  ! slot of 1024 'I' for next idle thread
5042  !$OMP PARALLEL DO &
5043  !$OMP PRIVATE(J,IJ) &
5044  !$OMP REDUCTION(+:B) &
5045  !$OMP SCHEDULE(DYNAMIC,ichunk)
5046  DO i=1,n
5047  ij=globalrowoffsets(i+l)+l
5048  b(i)=globalmatd(ij+i)*x(i)
5049  DO j=1,i-1
5050  b(j)=b(j)+globalmatd(ij+j)*x(i)
5051  b(i)=b(i)+globalmatd(ij+j)*x(j)
5052  END DO
5053  END DO
5054  !$OMP END PARALLEL DO
5055  ELSE
5056  ! sparse, compressed matrix
5057  IF(sparsematrixoffsets(2,1) /= n) THEN
5058  CALL peend(24,'Aborted, vector/matrix size mismatch')
5059  stop 'AVPRD0: mismatched vector and matrix'
5060  END IF
5061  ! parallelize row (group) loop
5062  ! slot of 1024 'I' for next idle thread
5063  !$OMP PARALLEL DO &
5064  !$OMP PRIVATE(I,IR,K,KK,LL,KU,INDID,INDIJ,J,JN,LJ) &
5065  !$OMP PRIVATE(IA,IB,IN,JA,JB) &
5066  !$OMP REDUCTION(+:B) &
5067  !$OMP SCHEDULE(DYNAMIC,ichunk)
5068  DO ipg=1,napgrp
5069  iproc=0
5070  !$ IPROC=OMP_GET_THREAD_NUM() ! thread number
5071  ! row group
5072  ia=globalallindexgroups(ipg) ! first (global) row
5073  ib=globalallindexgroups(ipg+1)-1 ! last (global) row
5074  in=ib-ia+1 ! number of rows
5075  !
5076  ! diagonal elements
5077  b(ia:ib)=globalmatd(ia:ib)*x(ia:ib)
5078  ! off-diagonals double precision
5079  ir=ipg
5080  kk=sparsematrixoffsets(1,ir) ! offset in 'd' (column lists)
5081  ll=sparsematrixoffsets(2,ir) ! offset in 'j' (matrix)
5082  ku=sparsematrixoffsets(1,ir+1)-kk
5083  indid=kk
5084  indij=ll
5085  IF (sparsematrixcolumns(indid+1) /= 0) THEN ! no compression
5086  DO i=ia,ib
5087  DO k=1,ku
5088  j=sparsematrixcolumns(indid+k)
5089  b(j)=b(j)+globalmatd(indij+k)*x(i)
5090  b(i)=b(i)+globalmatd(indij+k)*x(j)
5091  END DO
5092  indij=indij+ku
5093  END DO
5094  ELSE
5095  ! regions of continous column groups
5096  DO k=2,ku-2,2
5097  j=sparsematrixcolumns(indid+k) ! first group
5098  ja=globalallindexgroups(j) ! first (global) column
5099  lj=sparsematrixcolumns(indid+k-1) ! region offset
5100  jn=sparsematrixcolumns(indid+k+1)-lj ! number of columns
5101  jb=ja+jn-1 ! last (global) column
5102  lj=1 ! index (in group region)
5103  DO i=ia,ib
5104  b(i)=b(i)+dot_product(globalmatd(indij+lj:indij+lj+jn-1),x(ja:jb))
5105  lj=lj+jn
5106  END DO
5107  IF (mextnd == 0) THEN
5108  lj=1
5109  DO j=ja,jb
5110  b(j)=b(j)+dot_product(globalmatd(indij+lj:indij+jn*in:jn),x(ia:ib))
5111  lj=lj+1
5112  END DO
5113  END IF
5114  indij=indij+in*jn
5115  END DO
5116  END IF
5117  ! mixed precision
5118  IF (nspc > 1) THEN
5119  ir=ipg+napgrp+1 ! off-diagonals single precision
5120  kk=sparsematrixoffsets(1,ir) ! offset in 'd' (column lists)
5121  ll=sparsematrixoffsets(2,ir) ! offset in 'j' (matrix)
5122  ku=sparsematrixoffsets(1,ir+1)-kk
5123  indid=kk
5124  indij=ll
5125  IF (sparsematrixcolumns(indid+1) /= 0) THEN ! no compression
5126  DO i=ia,ib
5127  DO k=1,ku
5128  j=sparsematrixcolumns(indid+k)
5129  b(j)=b(j)+real(globalmatf(indij+k),mpd)*x(i)
5130  b(i)=b(i)+real(globalmatf(indij+k),mpd)*x(j)
5131  END DO
5132  indij=indij+ku
5133  END DO
5134  ELSE
5135  ! regions of continous column groups
5136  DO k=2,ku-2,2
5137  j=sparsematrixcolumns(indid+k) ! first group
5138  ja=globalallindexgroups(j) ! first (global) column
5139  lj=sparsematrixcolumns(indid+k-1) ! region offset
5140  jn=sparsematrixcolumns(indid+k+1)-lj ! number of columns
5141  jb=ja+jn-1 ! last (global) column
5142  lj=1 ! index (in group region)
5143  DO i=ia,ib
5144  b(i)=b(i)+dot_product(real(globalmatf(indij+lj:indij+lj+jn-1),mpd),x(ja:jb))
5145  lj=lj+jn
5146  END DO
5147  IF (mextnd == 0) THEN
5148  lj=1
5149  DO j=ja,jb
5150  b(j)=b(j)+dot_product(real(globalmatf(indij+lj:indij+jn*in:jn),mpd),x(ia:ib))
5151  lj=lj+1
5152  END DO
5153  END IF
5154  indij=indij+in*jn
5155  END DO
5156  END IF
5157  END IF
5158  END DO
5159  ENDIF
5160 
5161 END SUBROUTINE avprd0
5162 
5165 SUBROUTINE anasps
5166  USE mpmod
5167 
5168  IMPLICIT NONE
5169  INTEGER(mpi) :: ia
5170  INTEGER(mpi) :: ib
5171  INTEGER(mpi) :: in
5172  INTEGER(mpi) :: ipg
5173  INTEGER(mpi) :: ir
5174  INTEGER(mpi) :: ispc
5175  INTEGER(mpi) :: jn
5176  INTEGER(mpi) :: lj
5177  REAL(mps) :: avg
5178 
5179 
5180  INTEGER(mpl) :: k
5181  INTEGER(mpl) :: kk
5182  INTEGER(mpl) :: ku
5183  INTEGER(mpl) :: ll
5184  INTEGER(mpl) :: indid
5185  INTEGER(mpl), DIMENSION(12) :: icount
5186  SAVE
5187 
5188  ! require sparse storage
5189  IF(matsto /= 2) RETURN
5190  ! reset
5191  icount=0
5192  icount(4)=huge(icount(4))
5193  icount(7)=huge(icount(7))
5194  icount(10)=huge(icount(10))
5195  ! loop over precisions
5196  DO ispc=1,nspc
5197  ! loop over row groups
5198  DO ipg=1,napgrp
5199  ! row group
5200  ia=globalallindexgroups(ipg) ! first (global) row
5201  ib=globalallindexgroups(ipg+1)-1 ! last (global) row
5202  in=ib-ia+1 ! number of rows
5203 
5204  ir=ipg+(ispc-1)*(napgrp+1)
5205  kk=sparsematrixoffsets(1,ir) ! offset in 'd' (column lists)
5206  ll=sparsematrixoffsets(2,ir) ! offset in 'j' (matrix)
5207  ku=sparsematrixoffsets(1,ir+1)-kk
5208  indid=kk
5209  IF (sparsematrixcolumns(indid+1) /= 0) THEN ! no compression
5210  icount(1)=icount(1)+in
5211  icount(2)=icount(2)+in*ku
5212  ELSE
5213  ! regions of continous column groups
5214  DO k=2,ku-2,2
5215  lj=sparsematrixcolumns(indid+k-1) ! region offset
5216  jn=sparsematrixcolumns(indid+k+1)-lj ! number of columns
5217  icount(3)=icount(3)+1 ! block (region) counter
5218  icount(4)=min(icount(4),jn) ! min number of columns per block (region)
5219  icount(5)=icount(5)+jn ! sum number of columns per block (region)
5220  icount(6)=max(icount(6),jn) ! max number of columns per block (region)
5221  icount(7)=min(icount(7),in) ! min number of rows per block (region)
5222  icount(8)=icount(8)+in ! sum number of rows per block (region)
5223  icount(9)=max(icount(9),in) ! max number of rows per block (region)
5224  icount(10)=min(icount(10),in*jn) ! min number of elements per block (region)
5225  icount(11)=icount(11)+in*jn ! sum number of elements per block (region)
5226  icount(12)=max(icount(12),in*jn) ! max number of elements per block (region)
5227  END DO
5228  END IF
5229  END DO
5230  END DO
5231 
5232  WRITE(*,*) "analysis of sparsity structure"
5233  IF (icount(1) > 0) THEN
5234  WRITE(*,101) "rows without compression/blocks ", icount(1)
5235  WRITE(*,101) " contained elements ", icount(2)
5236  ENDIF
5237  WRITE(*,101) "number of block matrices ", icount(3)
5238  avg=real(icount(5),mps)/real(icount(3),mps)
5239  WRITE(*,101) "number of columns (min,mean,max) ", icount(4), avg, icount(6)
5240  avg=real(icount(8),mps)/real(icount(3),mps)
5241  WRITE(*,101) "number of rows (min,mean,max) ", icount(7), avg, icount(9)
5242  avg=real(icount(11),mps)/real(icount(3),mps)
5243  WRITE(*,101) "number of elements (min,mean,max) ", icount(10), avg, icount(12)
5244 101 FORMAT(2x,a34,i10,f10.3,i10)
5245 
5246 END SUBROUTINE anasps
5247 
5257 
5258 SUBROUTINE avprod(n,x,b)
5259  USE mpmod
5260 
5261  IMPLICIT NONE
5262 
5263  INTEGER(mpi), INTENT(IN) :: n
5264  REAL(mpd), INTENT(IN) :: x(n)
5265  REAL(mpd), INTENT(OUT) :: b(n)
5266 
5267  SAVE
5268  ! ...
5269  IF(n > nagb) THEN
5270  CALL peend(24,'Aborted, vector/matrix size mismatch')
5271  stop 'AVPROD: mismatched vector and matrix'
5272  END IF
5273  ! input to AVPRD0
5274  vecxav(1:n)=x
5275  vecxav(n+1:nagb)=0.0_mpd
5276  !use elimination for constraints ?
5277  IF(n < nagb) CALL qlmlq(vecxav,1,.false.) ! Q*x
5278  ! calclulate vecBav=globalMat*vecXav
5279  CALL avprd0(nagb,0_mpl,vecxav,vecbav)
5280  !use elimination for constraints ?
5281  IF(n < nagb) CALL qlmlq(vecbav,1,.true.) ! Q^t*x
5282  ! output from AVPRD0
5283  b=vecbav(1:n)
5284 
5285 END SUBROUTINE avprod
5286 
5287 
5297 
5298 SUBROUTINE ijpgrp(itema,itemb,ij,lr,iprc)
5299  USE mpmod
5300 
5301  IMPLICIT NONE
5302  INTEGER(mpi) :: ispc
5303  INTEGER(mpi) :: item1
5304  INTEGER(mpi) :: item2
5305  INTEGER(mpi) :: itemc
5306  INTEGER(mpi) :: jtem
5307  INTEGER(mpi) :: jtemn
5308  INTEGER(mpi) :: np
5309 
5310  INTEGER(mpi), INTENT(IN) :: itema
5311  INTEGER(mpi), INTENT(IN) :: itemb
5312  INTEGER(mpl), INTENT(OUT) :: ij
5313  INTEGER(mpi), INTENT(OUT) :: lr
5314  INTEGER(mpi), INTENT(OUT) :: iprc
5315 
5316  INTEGER(mpl) :: k
5317  INTEGER(mpl) :: kk
5318  INTEGER(mpl) :: kl
5319  INTEGER(mpl) :: ku
5320  INTEGER(mpl) :: ll
5321  ! ...
5322  ij=0
5323  lr=0
5324  iprc=0
5325  item1=max(itema,itemb) ! larger index
5326  item2=min(itema,itemb) ! smaller index
5327  IF(item2 <= 0.OR.item1 > napgrp) RETURN
5328  np=globalallindexgroups(item1+1)-globalallindexgroups(item1) ! size of group item1
5329  ! loop over precisions
5330  outer: DO ispc=1,nspc
5331  kk=sparsematrixoffsets(1,item1) ! offset (column lists)
5332  ll=sparsematrixoffsets(2,item1) ! offset (matrix)
5333  kl=1
5334  ku=sparsematrixoffsets(1,item1+1)-kk
5335  item1=item1+napgrp+1
5336  iprc=ispc
5337  IF (sparsematrixcolumns(kk+1) == 0) THEN ! compression ?
5338  ! compressed (list of continous regions of parameter groups (pairs of offset and 1. group index)
5339  kl=2
5340  ku=ku-2
5341  IF(ku < kl) cycle outer ! not found
5342  DO
5343  k=2*((kl+ku)/4) ! binary search
5344  jtem=sparsematrixcolumns(kk+k) ! first column (group) of region
5345  jtemn=sparsematrixcolumns(kk+k+2) ! first column (group) after region
5346  IF(item2 >= jtem.AND.item2 < jtemn) THEN
5347  ! length of region
5348  lr=sparsematrixcolumns(kk+k+1)-sparsematrixcolumns(kk+k-1)
5349  IF (globalallindexgroups(item2)-globalallindexgroups(jtem) >= lr) cycle outer ! outside region
5350  EXIT ! found
5351  END IF
5352  IF(item2 < jtem) THEN
5353  ku=k-2
5354  ELSE IF(item2 >= jtemn) THEN
5355  kl=k+2
5356  END IF
5357  IF(kl <= ku) cycle
5358  cycle outer ! not found
5359  END DO
5360  ! group offset in row
5361  ij=sparsematrixcolumns(kk+k-1)
5362  ! absolute offset
5363  ij=ll+ij*np+globalallindexgroups(item2)-globalallindexgroups(jtem)+1
5364 
5365  ELSE
5366  ! simple column list
5367  itemc=globalallindexgroups(item2) ! first (col) index of group
5368  lr=int(ku,mpi) ! number of columns
5369  IF(ku < kl) cycle outer ! not found
5370  DO
5371  k=(kl+ku)/2 ! binary search
5372  jtem=sparsematrixcolumns(kk+k)
5373  IF(itemc == jtem) EXIT ! found
5374  IF(itemc < jtem) THEN
5375  ku=k-1
5376  ELSE IF(itemc > jtem) THEN
5377  kl=k+1
5378  END IF
5379  IF(kl <= ku) cycle
5380  cycle outer ! not found
5381  END DO
5382  ij=ll+k
5383 
5384  END IF
5385  RETURN
5386  END DO outer
5387 
5388 END SUBROUTINE ijpgrp
5389 
5395 
5396 FUNCTION ijprec(itema,itemb)
5397  USE mpmod
5398 
5399  IMPLICIT NONE
5400 
5401  INTEGER(mpi) :: lr
5402  INTEGER(mpl) :: ij
5403 
5404  INTEGER(mpi), INTENT(IN) :: itema
5405  INTEGER(mpi), INTENT(IN) :: itemb
5406  INTEGER(mpi) :: ijprec
5407 
5408  ! ...
5409  ijprec=1
5410  IF (matsto == 2.AND.nspc > 1) THEN ! sparse storage with mixed precision
5411  ! check groups
5412  CALL ijpgrp(itema,itemb,ij,lr,ijprec)
5413  END IF
5414 
5415 END FUNCTION ijprec
5416 
5424 
5425 FUNCTION ijadd(itema,itemb) ! index using "d" and "z"
5426  USE mpmod
5427 
5428  IMPLICIT NONE
5429 
5430  INTEGER(mpi) :: item1
5431  INTEGER(mpi) :: item2
5432  INTEGER(mpi) :: ipg1
5433  INTEGER(mpi) :: ipg2
5434  INTEGER(mpi) :: lr
5435  INTEGER(mpi) :: iprc
5436 
5437  INTEGER(mpi), INTENT(IN) :: itema
5438  INTEGER(mpi), INTENT(IN) :: itemb
5439 
5440  INTEGER(mpl) :: ijadd
5441  INTEGER(mpl) :: ij
5442  ! ...
5443  ijadd=0
5444  item1=max(itema,itemb) ! larger index
5445  item2=min(itema,itemb) ! smaller index
5446  !print *, ' ijadd ', item1, item2
5447  IF(item2 <= 0.OR.item1 > nagb) RETURN
5448  IF(item1 == item2) THEN ! diagonal element
5449  ijadd=item1
5450  RETURN
5451  END IF
5452  ! ! off-diagonal element
5453  ! get parameter groups
5454  ipg1=globalallpartogroup(item1)
5455  ipg2=globalallpartogroup(item2)
5456  ! get offset for groups
5457  CALL ijpgrp(ipg1,ipg2,ij,lr,iprc)
5458  IF (ij == 0) RETURN
5459  ! add offset inside groups
5460  ijadd=ij+(item2-globalallindexgroups(ipg2))+(item1-globalallindexgroups(ipg1))*lr
5461  ! reduced precision?
5462  IF (iprc > 1) ijadd=-ijadd
5463 
5464 END FUNCTION ijadd
5465 
5471 
5472 FUNCTION matij(itema,itemb)
5473  USE mpmod
5474 
5475  IMPLICIT NONE
5476 
5477  INTEGER(mpi) :: item1
5478  INTEGER(mpi) :: item2
5479  INTEGER(mpl) :: i
5480  INTEGER(mpl) :: j
5481  INTEGER(mpl) :: ij
5482  INTEGER(mpl) :: ijadd
5483 
5484  INTEGER(mpi), INTENT(IN) :: itema
5485  INTEGER(mpi), INTENT(IN) :: itemb
5486 
5487  REAL(mpd) :: matij
5488  ! ...
5489  matij=0.0_mpd
5490  item1=max(itema,itemb) ! larger index
5491  item2=min(itema,itemb) ! smaller index
5492  IF(item2 <= 0.OR.item1 > nagb) RETURN
5493 
5494  i=item1
5495  j=item2
5496 
5497  IF(matsto /= 2) THEN ! full or unpacked (block diagonal) symmetric matrix
5498  ij=globalrowoffsets(i)+j
5499  matij=globalmatd(ij)
5500  ELSE ! sparse symmetric matrix
5501  ij=ijadd(item1,item2) ! inline code requires same time
5502  IF (ij > 0) THEN
5503  matij=globalmatd(ij)
5504  ELSE IF (ij < 0) THEN
5505  matij=real(globalmatf(-ij),mpd)
5506  END IF
5507  END IF
5508 
5509 END FUNCTION matij
5510 
5513 
5514 SUBROUTINE mhalf2
5515  USE mpmod
5516 
5517  IMPLICIT NONE
5518  INTEGER(mpi) :: i
5519  INTEGER(mpi) :: ia
5520  INTEGER(mpi) :: ib
5521  INTEGER(mpi) :: ichunk
5522  INTEGER(mpi) :: in
5523  INTEGER(mpi) :: ipg
5524  INTEGER(mpi) :: ir
5525  INTEGER(mpi) :: ispc
5526  INTEGER(mpi) :: j
5527  INTEGER(mpi) :: ja
5528  INTEGER(mpi) :: jb
5529  INTEGER(mpi) :: jn
5530  INTEGER(mpi) :: lj
5531 
5532  INTEGER(mpl) :: ij
5533  INTEGER(mpl) :: ijadd
5534  INTEGER(mpl) :: k
5535  INTEGER(mpl) :: kk
5536  INTEGER(mpl) :: ku
5537  INTEGER(mpl) :: ll
5538  ! ...
5539 
5540  ichunk=min((napgrp+mthrd-1)/mthrd/8+1,1024)
5541 
5542  DO ispc=1,nspc
5543  ! parallelize row loop
5544  ! slot of 1024 'I' for next idle thread
5545  !$OMP PARALLEL DO &
5546  !$OMP PRIVATE(I,IR,K,KK,LL,KU,IJ,J,LJ) &
5547  !$OMP PRIVATE(IA,IB,IN,JA,JB,JN) &
5548  !$OMP SCHEDULE(DYNAMIC,ichunk)
5549  DO ipg=1,napgrp
5550  ! row group
5551  ia=globalallindexgroups(ipg) ! first (global) row
5552  ib=globalallindexgroups(ipg+1)-1 ! last (global) row
5553  in=ib-ia+1 ! number of rows
5554  !
5555  ir=ipg+(ispc-1)*(napgrp+1)
5556  kk=sparsematrixoffsets(1,ir) ! offset in 'd' (column lists)
5557  ll=sparsematrixoffsets(2,ir) ! offset in 'j' (matrix)
5558  ku=sparsematrixoffsets(1,ir+1)-kk
5559  ! regions of continous column groups
5560  DO k=2,ku-2,2
5561  j=sparsematrixcolumns(kk+k) ! first group
5562  ja=globalallindexgroups(j) ! first (global) column
5563  lj=sparsematrixcolumns(kk+k-1) ! region offset
5564  jn=sparsematrixcolumns(kk+k+1)-lj ! number of columns
5565  jb=ja+jn-1 ! last (global) column
5566  ! skip first half
5567  IF (sparsematrixcolumns(kk+k+2) <= ipg) THEN
5568  ll=ll+in*jn
5569  cycle
5570  END IF
5571  ! at diagonal or in second half
5572  DO i=ia,ib ! loop over rows
5573  DO j=ja,jb ! loop over columns
5574  ll=ll+1
5575  IF (j > i) THEN
5576  ij=ijadd(i,j)
5577  IF (ispc==1) THEN
5578  globalmatd(ll)=globalmatd(ij)
5579  ELSE
5580  globalmatf(ll)=globalmatf(-ij)
5581  END IF
5582  END IF
5583  END DO
5584  END DO
5585  END DO
5586  END DO
5587  !$OMP END PARALLEL DO
5588  END DO
5589 
5590 END SUBROUTINE mhalf2
5591 
5600 
5601 SUBROUTINE sechms(deltat,nhour,minut,secnd)
5602  USE mpdef
5603 
5604  IMPLICIT NONE
5605  REAL(mps), INTENT(IN) :: deltat
5606  INTEGER(mpi), INTENT(OUT) :: minut
5607  INTEGER(mpi), INTENT(OUT):: nhour
5608  REAL(mps), INTENT(OUT):: secnd
5609  INTEGER(mpi) :: nsecd
5610  ! DELTAT = time in sec -> NHOUR,MINUT,SECND
5611  ! ...
5612  nsecd=nint(deltat,mpi) ! -> integer
5613  nhour=nsecd/3600
5614  minut=nsecd/60-60*nhour
5615  secnd=deltat-60*(minut+60*nhour)
5616 END SUBROUTINE sechms
5617 
5645 
5646 INTEGER(mpi) FUNCTION inone(item) ! translate 1-D identifier to nrs
5647  USE mpmod
5648  USE mpdalc
5649 
5650  IMPLICIT NONE
5651  INTEGER(mpi), INTENT(IN) :: item
5652  INTEGER(mpi) :: j
5653  INTEGER(mpi) :: k
5654  INTEGER(mpi) :: iprime
5655  INTEGER(mpl) :: length
5656  INTEGER(mpl), PARAMETER :: four = 4
5657 
5658  inone=0
5659  !print *, ' INONE ', item
5660  IF(item <= 0) RETURN
5661  IF(globalparheader(-1) == 0) THEN
5662  length=128 ! initial number
5663  CALL mpalloc(globalparlabelindex,four,length,'INONE: label & index')
5664  CALL mpalloc(globalparhashtable,2*length,'INONE: hash pointer')
5665  globalparhashtable = 0
5666  globalparheader(-0)=int(length,mpi) ! length of labels/indices
5667  globalparheader(-1)=0 ! number of stored items
5668  globalparheader(-2)=0 ! =0 during build-up
5669  globalparheader(-3)=int(length,mpi) ! next number
5670  globalparheader(-4)=iprime(globalparheader(-0)) ! prime number
5671  globalparheader(-5)=0 ! number of overflows
5672  globalparheader(-6)=0 ! nr of variable parameters
5673  globalparheader(-8)=0 ! number of sorted items
5674  END IF
5675  outer: DO
5676  j=1+mod(item,globalparheader(-4))+globalparheader(-0)
5677  inner: DO ! normal case: find item
5678  k=j
5679  j=globalparhashtable(k)
5680  IF(j == 0) EXIT inner ! unused hash code
5681  IF(item == globalparlabelindex(1,j)) EXIT outer ! found
5682  END DO inner
5683  ! not found
5684  IF(globalparheader(-1) == globalparheader(-0).OR.globalparheader(-2) /= 0) THEN
5685  globalparheader(-5)=globalparheader(-5)+1 ! overflow
5686  j=0
5687  RETURN
5688  END IF
5689  globalparheader(-1)=globalparheader(-1)+1 ! increase number of elements
5691  j=globalparheader(-1)
5692  globalparhashtable(k)=j ! hash index
5693  globalparlabelindex(1,j)=item ! add new item
5694  globalparlabelindex(2,j)=0 ! reset counter
5695  globalparlabelindex(3,j)=0 ! reset group info
5696  globalparlabelindex(4,j)=0 ! reset group info
5697  IF(globalparheader(-1) /= globalparheader(-0)) EXIT outer
5698  ! update with larger dimension and redefine index
5700  CALL upone
5701  IF (lvllog > 1) WRITE(lunlog,*) 'INONE: array increased to', &
5702  globalparheader(-3),' words'
5703  END DO outer
5704 
5705  IF(globalparheader(-2) == 0) THEN
5706  globalparlabelindex(2,j)=globalparlabelindex(2,j)+1 ! increase counter
5708  END IF
5709  inone=j
5710 END FUNCTION inone
5711 
5713 SUBROUTINE upone
5714  USE mpmod
5715  USE mpdalc
5716 
5717  IMPLICIT NONE
5718  INTEGER(mpi) :: i
5719  INTEGER(mpi) :: j
5720  INTEGER(mpi) :: k
5721  INTEGER(mpi) :: iprime
5722  INTEGER(mpi) :: nused
5723  LOGICAL :: finalUpdate
5724  INTEGER(mpl) :: oldLength
5725  INTEGER(mpl) :: newLength
5726  INTEGER(mpl), PARAMETER :: four = 4
5727  INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: tempArr
5728  SAVE
5729  ! ...
5730  finalupdate=(globalparheader(-3) == globalparheader(-1))
5731  IF(finalupdate) THEN ! final (cleanup) call
5732  IF (globalparheader(-1) > globalparheader(-8)) THEN
5733  CALL sort22(globalparlabelindex,globalparheader(-1)) ! sort items
5735  END IF
5736  END IF
5737  ! save old LabelIndex
5738  nused = globalparheader(-1)
5739  oldlength = globalparheader(-0)
5740  CALL mpalloc(temparr,four,oldlength,'INONE: temp array')
5741  temparr(:,1:nused)=globalparlabelindex(:,1:nused)
5744  ! create new LabelIndex
5745  newlength = globalparheader(-3)
5746  CALL mpalloc(globalparlabelindex,four,newlength,'INONE: label & index')
5747  CALL mpalloc(globalparhashtable,2*newlength,'INONE: hash pointer')
5748  globalparhashtable = 0
5749  globalparlabelindex(:,1:nused) = temparr(:,1:nused) ! copy back saved content
5750  CALL mpdealloc(temparr)
5751  globalparheader(-0)=int(newlength,mpi) ! length of labels/indices
5753  globalparheader(-4)=iprime(globalparheader(-0)) ! prime number < LNDA
5754  ! redefine hash
5755  outer: DO i=1,globalparheader(-1)
5757  inner: DO
5758  k=j
5759  j=globalparhashtable(k)
5760  IF(j == 0) EXIT inner ! unused hash code
5761  IF(j == i) cycle outer ! found
5762  ENDDO inner
5763  globalparhashtable(k)=i
5764  END DO outer
5765  IF(.NOT.finalupdate) RETURN
5766 
5767  globalparheader(-2)=1 ! set flag to inhibit further updates
5768  IF (lvllog > 1) THEN
5769  WRITE(lunlog,*) ' '
5770  WRITE(lunlog,*) 'INONE: array reduced to',newlength,' words'
5771  WRITE(lunlog,*) 'INONE:',globalparheader(-1),' items stored.'
5772  END IF
5773 END SUBROUTINE upone ! update, redefine
5774 
5776 SUBROUTINE useone
5777  USE mpmod
5778 
5779  IMPLICIT NONE
5780  INTEGER(mpi) :: i
5781  INTEGER(mpi) :: j
5782  INTEGER(mpi) :: k
5783  SAVE
5784  ! ...
5785  IF (globalparheader(-1) > globalparheader(-8)) THEN
5786  CALL sort22(globalparlabelindex,globalparheader(-1)) ! sort items
5787  ! redefine hash
5788  globalparhashtable = 0
5789  outer: DO i=1,globalparheader(-1)
5791  inner: DO
5792  k=j
5793  j=globalparhashtable(k)
5794  IF(j == 0) EXIT inner ! unused hash code
5795  IF(j == i) cycle outer ! found
5796  ENDDO inner
5797  globalparhashtable(k)=i
5798  END DO outer
5800  END IF
5801 END SUBROUTINE useone ! make usable
5802 
5807 
5808 INTEGER(mpi) FUNCTION iprime(n)
5809  USE mpdef
5810 
5811  IMPLICIT NONE
5812  INTEGER(mpi), INTENT(IN) :: n
5813  INTEGER(mpi) :: nprime
5814  INTEGER(mpi) :: nsqrt
5815  INTEGER(mpi) :: i
5816  ! ...
5817  SAVE
5818  nprime=n ! max number
5819  IF(mod(nprime,2) == 0) nprime=nprime+1 ! ... odd number
5820  outer: DO
5821  nprime=nprime-2 ! next lower odd number
5822  nsqrt=int(sqrt(real(nprime,mps)),mpi)
5823  DO i=3,nsqrt,2 !
5824  IF(i*(nprime/i) == nprime) cycle outer ! test prime number
5825  END DO
5826  EXIT outer ! found
5827  END DO outer
5828  iprime=nprime
5829 END FUNCTION iprime
5830 
5840 SUBROUTINE loop1
5841  USE mpmod
5842  USE mpdalc
5843 
5844  IMPLICIT NONE
5845  INTEGER(mpi) :: i
5846  INTEGER(mpi) :: idum
5847  INTEGER(mpi) :: in
5848  INTEGER(mpi) :: indab
5849  INTEGER(mpi) :: itgbi
5850  INTEGER(mpi) :: itgbl
5851  INTEGER(mpi) :: ivgbi
5852  INTEGER(mpi) :: j
5853  INTEGER(mpi) :: jgrp
5854  INTEGER(mpi) :: lgrp
5855  INTEGER(mpi) :: mqi
5856  INTEGER(mpi) :: nc31
5857  INTEGER(mpi) :: nr
5858  INTEGER(mpi) :: nwrd
5859  INTEGER(mpi) :: inone
5860  REAL(mpd) :: param
5861  REAL(mpd) :: presg
5862  REAL(mpd) :: prewt
5863 
5864  INTEGER(mpl) :: length
5865  INTEGER(mpl) :: rows
5866  SAVE
5867  ! ...
5868  WRITE(lunlog,*) ' '
5869  WRITE(lunlog,*) 'LOOP1: starting'
5870  CALL mstart('LOOP1')
5871 
5872  ! add labels from parameter, constraints, measurements -------------
5873  DO i=1, lenparameters
5874  idum=inone(listparameters(i)%label)
5875  END DO
5876  DO i=1, lenpresigmas
5877  idum=inone(listpresigmas(i)%label)
5878  END DO
5879  DO i=1, lenconstraints
5880  idum=inone(listconstraints(i)%label)
5881  END DO
5882  DO i=1, lenmeasurements
5883  idum=inone(listmeasurements(i)%label)
5884  END DO
5885 
5886  IF(globalparheader(-1) /= 0) THEN
5887  WRITE(lunlog,*) 'LOOP1:',globalparheader(-1), ' labels from txt data stored'
5888  END IF
5889  WRITE(lunlog,*) 'LOOP1: reading data files'
5890 
5891  DO
5892  DO j=1,globalparheader(-1)
5893  globalparlabelindex(2,j)=0 ! reset count
5894  END DO
5895 
5896  ! read all data files and add all labels to global labels table ----
5897 
5898  IF(mprint /= 0) THEN
5899  WRITE(*,*) 'Read all binary data files:'
5900  END IF
5901  CALL hmpldf(1,'Number of words/record in binary file')
5902  CALL hmpdef(8,0.0,60.0,'not_stored data per record')
5903  ! define read buffer
5904  nc31=ncache/(31*mthrdr) ! split read cache 1 : 10 : 10*2 for pointers, ints, floats
5905  nwrd=nc31+1
5906  length=nwrd*mthrdr
5907  CALL mpalloc(readbufferpointer,length,'read buffer, pointer')
5908  nwrd=nc31*10+2+ndimbuf
5909  length=nwrd*mthrdr
5910  CALL mpalloc(readbufferdatai,length,'read buffer, integer')
5911  CALL mpalloc(readbufferdatad,length,'read buffer, double')
5912  ! to read (old) float binary files
5913  length=(ndimbuf+2)*mthrdr
5914  CALL mpalloc(readbufferdataf,length,'read buffer, float')
5915 
5916  DO
5917  CALL peread(nr) ! read records
5918  IF (skippedrecords == 0) THEN
5919  CALL peprep(0) ! prepare records
5920  CALL pepgrp ! update parameter group info
5921  END IF
5922  IF(nr <= 0) EXIT ! end of data?
5923  END DO
5924  ! release read buffer
5929  IF (skippedrecords == 0) THEN
5930  EXIT
5931  ELSE
5932  WRITE(lunlog,*) 'LOOP1: reading data files again'
5933  END IF
5934  END DO
5935 
5936  IF(nhistp /= 0) THEN
5937  CALL hmprnt(1)
5938  CALL hmprnt(8)
5939  END IF
5940  CALL hmpwrt(1)
5941  CALL hmpwrt(8)
5942  ntgb = globalparheader(-1) ! total number of labels/parameters
5943  IF (ntgb == 0) THEN
5944  CALL peend(21,'Aborted, no labels/parameters defined')
5945  stop 'LOOP1: no labels/parameters defined'
5946  END IF
5947  CALL upone ! finalize the global label table
5948 
5949  WRITE(lunlog,*) 'LOOP1:',ntgb, &
5950  ' is total number NTGB of labels/parameters'
5951  ! histogram number of entries per label ----------------------------
5952  CALL hmpldf(2,'Number of entries per label')
5953  DO j=1,ntgb
5954  CALL hmplnt(2,globalparlabelindex(2,j))
5955  END DO
5956  IF(nhistp /= 0) CALL hmprnt(2) ! print histogram
5957  CALL hmpwrt(2) ! write to his file
5958 
5959  ! three subarrays for all global parameters ------------------------
5960  length=ntgb
5961  CALL mpalloc(globalparameter,length,'global parameters')
5962  globalparameter=0.0_mpd
5963  CALL mpalloc(globalparpresigma,length,'pre-sigmas') ! presigmas
5965  CALL mpalloc(globalparstart,length,'global parameters at start')
5966  globalparstart=0.
5967  CALL mpalloc(globalparcopy,length,'copy of global parameters')
5968  CALL mpalloc(globalparcounts,length,'global parameter counts')
5969  CALL mpalloc(globalparcons,length,'global parameter constraints')
5970  globalparcons=0
5971 
5972  DO i=1,lenparameters ! parameter start values
5973  param=listparameters(i)%value
5974  in=inone(listparameters(i)%label)
5975  IF(in /= 0) THEN
5976  globalparameter(in)=param
5977  globalparstart(in)=param
5978  ENDIF
5979  END DO
5980 
5981  npresg=0
5982  DO i=1,lenpresigmas ! pre-sigma values
5983  presg=listpresigmas(i)%value
5984  in=inone(listpresigmas(i)%label)
5985  IF(in /= 0) THEN
5986  IF(presg > 0.0) npresg=npresg+1 ! FIXME: check if enough 'entries'?
5987  globalparpresigma(in)=presg ! insert pre-sigma 0 or > 0
5988  END IF
5989  END DO
5990  WRITE(lunlog,*) 'LOOP1:',npresg,' is number of pre-sigmas'
5991  WRITE(*,*) 'LOOP1:',npresg,' is number of pre-sigmas'
5992  IF(npresg == 0) WRITE(*,*) 'Warning: no pre-sigmas defined'
5993 
5994  ! determine flag variable (active) or fixed (inactive) -------------
5995 
5996  indab=0
5997  DO i=1,ntgb
5999  IF (globalparpresigma(i) < 0.0) THEN
6000  globalparlabelindex(2,i)=-1 ! fixed (pre-sigma), not used in matrix (not active)
6001  ELSE IF(globalparcounts(i) < mreqenf) THEN
6002  globalparlabelindex(2,i)=-2 ! fixed (entries cut), not used in matrix (not active)
6003  ELSE
6004  indab=indab+1
6005  globalparlabelindex(2,i)=indab ! variable, used in matrix (active)
6006  END IF
6007  END DO
6008  globalparheader(-6)=indab ! counted variable
6009  nvgb=indab ! nr of variable parameters
6010  WRITE(lunlog,*) 'LOOP1:',nvgb, ' is number NVGB of variable parameters'
6011  IF(iteren > mreqenf) THEN
6012  IF (mcount == 0) THEN
6013  CALL loop1i ! iterate entries cut
6014  ELSE
6015  WRITE(lunlog,*) 'LOOP1: counting records, NO iteration of entries cut !'
6016  iteren=0
6017  END IF
6018  END IF
6019 
6020  ! --- check for parameter groups
6021  CALL hmpdef(15,0.0,120.0,'Number of parameters per group')
6022  ntpgrp=0
6023  DO j=1,ntgb
6024  IF (globalparlabelindex(3,j) == 0) cycle ! skip empty parameter
6025  ! new group?
6027  globalparlabelindex(4,j)=ntpgrp ! relation total index -> group
6028  END DO
6029  ! check variable parameters
6030  nvpgrp=0
6031  lgrp=-1
6032  DO j=1,ntgb
6033  IF (globalparlabelindex(2,j) <= 0) cycle ! skip fixed parameter
6034  ! new group ?
6035  IF (globalparlabelindex(4,j) /= lgrp) nvpgrp=nvpgrp+1
6036  lgrp=globalparlabelindex(4,j)
6037  END DO
6038  length=ntpgrp; rows=2
6039  CALL mpalloc(globaltotindexgroups,rows,length,'parameter groups, 1. index and size')
6041  ! fill
6042  lgrp=-1
6043  DO j=1,ntgb
6044  IF (globalparlabelindex(3,j) == 0) cycle ! skip empty parameter
6045  jgrp=globalparlabelindex(4,j)
6046  IF (jgrp /= lgrp) globaltotindexgroups(1,jgrp)=j ! first (total) index
6047  globaltotindexgroups(2,jgrp)=globaltotindexgroups(2,jgrp)+1 ! (total) size
6048  lgrp=jgrp
6049  END DO
6050  DO j=1,ntpgrp
6051  CALL hmpent(15,real(globaltotindexgroups(2,j),mps))
6052  END DO
6053  IF(nhistp /= 0) CALL hmprnt(15) ! print histogram
6054  CALL hmpwrt(15) ! write to his file
6055  WRITE(lunlog,*) 'LOOP1:',ntpgrp, &
6056  ' is total number NTPGRP of label/parameter groups'
6057  !print *, ' globalTotIndexGroups ', globalTotIndexGroups
6058 
6059  ! translation table of length NVGB of total global indices ---------
6060  length=nvgb
6061  CALL mpalloc(globalparvartototal,length,'translation table var -> total')
6062  indab=0
6063  DO i=1,ntgb
6064  IF(globalparlabelindex(2,i) > 0) THEN
6065  indab=indab+1
6066  globalparvartototal(indab)=i
6067  END IF
6068  END DO
6069 
6070  ! regularization ---------------------------------------------------
6071  CALL mpalloc(globalparpreweight,length,'pre-sigmas weights') ! presigma weights
6072  WRITE(*,112) ' Default pre-sigma =',regpre, &
6073  ' (if no individual pre-sigma defined)'
6074  WRITE(*,*) 'Pre-sigma factor is',regula
6075 
6076  IF(nregul == 0) THEN
6077  WRITE(*,*) 'No regularization will be done'
6078  ELSE
6079  WRITE(*,*) 'Regularization will be done, using factor',regula
6080  END IF
6081 112 FORMAT(a,e9.2,a)
6082  IF (nvgb <= 0) THEN
6083  CALL peend(22,'Aborted, no variable global parameters')
6084  stop '... no variable global parameters'
6085  ENDIF
6086 
6087  DO ivgbi=1,nvgb ! IVGBI = variable parameter index
6088  itgbi=globalparvartototal(ivgbi) ! ITGBI = global parameter index
6089  presg=globalparpresigma(itgbi) ! get pre-sigma
6090  prewt=0.0 ! pre-weight
6091  IF(presg > 0.0) THEN
6092  prewt=1.0/presg**2 ! 1/presigma^2
6093  ELSE IF(presg == 0.0.AND.regpre > 0.0) THEN
6094  prewt=1.0/real(regpre**2,mpd) ! default 1/presigma^2
6095  END IF
6096  globalparpreweight(ivgbi)=regula*prewt ! weight = factor / presigma^2
6097  END DO
6098 
6099  ! WRITE(*,*) 'GlPa_index GlPa_label array1 array6'
6100  DO i=1,ntgb
6101  itgbl=globalparlabelindex(1,i)
6102  ivgbi=globalparlabelindex(2,i)
6103  IF(ivgbi > 0) THEN
6104  ! WRITE(*,111) I,ITGBL,QM(IND1+I),QM(IND6+IVGBI)
6105  ELSE
6106  ! WRITE(*,111) I,ITGBL,QM(IND1+I)
6107  END IF
6108  END DO
6109  ! 111 FORMAT(I5,I10,F10.5,E12.4)
6110  WRITE(*,101) 'NTGB',ntgb,'total number of parameters'
6111  WRITE(*,101) 'NVGB',nvgb,'number of variable parameters'
6112  ! To avoid INT(mpi) overflows in diagonalization
6113  IF (metsol == 2.AND.nvgb >= 46340) THEN
6114  metsol=1
6115  WRITE(*,101) 'Too many variable parameters for diagonalization, fallback is inversion'
6116  END IF
6117 
6118  ! print overview over important numbers ----------------------------
6119 
6120  nrecal=nrec
6121  IF(mprint /= 0) THEN
6122  WRITE(*,*) ' '
6123  WRITE(*,101) ' NREC',nrec,'number of records'
6124  IF (nrecd > 0) WRITE(*,101) ' NRECD',nrec,'number of records containing doubles'
6125  IF (mcount == 0) THEN
6126  WRITE(*,101) 'MREQENF',mreqenf,'required number of entries (eqns in binary files)'
6127  ELSE
6128  WRITE(*,101) 'MREQENF',mreqenf,'required number of entries (recs in binary files)'
6129  ENDIF
6130  IF(iteren > mreqenf) &
6131  WRITE(*,101) 'ITEREN',iteren,'iterate cut for parameters with less entries'
6132  WRITE(*,101) 'MREQENA',mreqena,'required number of entries (from accepted fits)'
6133  IF (mreqpe > 1) WRITE(*,101) &
6134  'MREQPE',mreqpe,'required number of pair entries'
6135  IF (msngpe >= 1) WRITE(*,101) &
6136  'MSNGPE',msngpe,'max pair entries single prec. storage'
6137  WRITE(*,101) 'NTGB',ntgb,'total number of parameters'
6138  WRITE(*,101) 'NVGB',nvgb,'number of variable parameters'
6139  IF(mprint > 1) THEN
6140  WRITE(*,*) ' '
6141  WRITE(*,*) 'Global parameter labels:'
6142  mqi=ntgb
6143  IF(mqi <= 100) THEN
6144  WRITE(*,*) (globalparlabelindex(2,i),i=1,mqi)
6145  ELSE
6146  WRITE(*,*) (globalparlabelindex(2,i),i=1,30)
6147  WRITE(*,*) ' ...'
6148  mqi=((mqi-20)/20)*20+1
6149  WRITE(*,*) (globalparlabelindex(2,i),i=mqi,ntgb)
6150  END IF
6151  END IF
6152  WRITE(*,*) ' '
6153  WRITE(*,*) ' '
6154  END IF
6155  WRITE(8,*) ' '
6156  WRITE(8,101) ' NREC',nrec,'number of records'
6157  IF (nrecd > 0) WRITE(8,101) ' NRECD',nrec,'number of records containing doubles'
6158  IF (mcount == 0) THEN
6159  WRITE(8,101) 'MREQENF',mreqenf,'required number of entries (eqns in binary files)'
6160  ELSE
6161  WRITE(8,101) 'MREQENF',mreqenf,'required number of entries (recs in binary files)'
6162  ENDIF
6163  IF(iteren > mreqenf) &
6164  WRITE(8,101) 'ITEREN',iteren,'iterate cut for parameters with less entries'
6165  WRITE(8,101) 'MREQENA',mreqena,'required number of entries (from accepted fits)'
6166 
6167  WRITE(lunlog,*) 'LOOP1: ending'
6168  WRITE(lunlog,*) ' '
6169  CALL mend
6170 
6171 101 FORMAT(1x,a8,' =',i10,' = ',a)
6172 END SUBROUTINE loop1
6173 
6181 SUBROUTINE loop1i
6182  USE mpmod
6183  USE mpdalc
6184 
6185  IMPLICIT NONE
6186  INTEGER(mpi) :: i
6187  INTEGER(mpi) :: ibuf
6188  INTEGER(mpi) :: ij
6189  INTEGER(mpi) :: indab
6190  INTEGER(mpi) :: inder
6191  INTEGER(mpi) :: isfrst
6192  INTEGER(mpi) :: islast
6193  INTEGER(mpi) :: ist
6194  INTEGER(mpi) :: j
6195  INTEGER(mpi) :: ja
6196  INTEGER(mpi) :: jb
6197  INTEGER(mpi) :: jsp
6198  INTEGER(mpi) :: nc31
6199  INTEGER(mpi) :: nr
6200  INTEGER(mpi) :: nlow
6201  INTEGER(mpi) :: nst
6202  INTEGER(mpi) :: nwrd
6203  REAL(mpr8) :: glder
6204 
6205  INTEGER(mpl) :: length
6206  INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: newCounter
6207  SAVE
6208 
6209  isfrst(ibuf)=readbufferpointer(ibuf)+1
6210  islast(ibuf)=readbufferdatai(readbufferpointer(ibuf))
6211  inder(i)=readbufferdatai(i)
6212  glder(i)=readbufferdatad(i)
6213  ! ...
6214  WRITE(lunlog,*) ' '
6215  WRITE(lunlog,*) 'LOOP1: iterating'
6216  WRITE(*,*) ' '
6217  WRITE(*,*) 'LOOP1: iterating'
6218 
6219  length=ntgb
6220  CALL mpalloc(newcounter,length,'new entries counter')
6221  newcounter=0
6222 
6223  ! define read buffer
6224  nc31=ncache/(31*mthrdr) ! split read cache 1 : 10 : 10*2 for pointers, ints, floats
6225  nwrd=nc31+1
6226  length=nwrd*mthrdr
6227  CALL mpalloc(readbufferpointer,length,'read buffer, pointer')
6228  nwrd=nc31*10+2+ndimbuf
6229  length=nwrd*mthrdr
6230  CALL mpalloc(readbufferdatai,length,'read buffer, integer')
6231  CALL mpalloc(readbufferdatad,length,'read buffer, double')
6232  ! to read (old) float binary files
6233  length=(ndimbuf+2)*mthrdr
6234  CALL mpalloc(readbufferdataf,length,'read buffer, float')
6235 
6236  DO
6237  CALL peread(nr) ! read records
6238  CALL peprep(1) ! prepare records
6239  DO ibuf=1,numreadbuffer ! buffer for current record
6240  ist=isfrst(ibuf)
6241  nst=islast(ibuf)
6242  nwrd=nst-ist+1
6243  DO ! loop over measurements
6244  CALL isjajb(nst,ist,ja,jb,jsp)
6245  IF(ja == 0.AND.jb == 0) EXIT
6246  IF(ja /= 0) THEN
6247  nlow=0
6248  DO j=1,ist-jb
6249  ij=inder(jb+j) ! index of global parameter
6250  ij=globalparlabelindex(2,ij) ! change to variable parameter
6251  IF(ij == -2) nlow=nlow+1 ! fixed by entries cut
6252  END DO
6253  IF(nlow == 0) THEN
6254  DO j=1,ist-jb
6255  ij=inder(jb+j) ! index of global parameter
6256  newcounter(ij)=newcounter(ij)+1 ! count again
6257  END DO
6258  ENDIF
6259  END IF
6260  END DO
6261  ! end-of-event
6262  END DO
6263  IF(nr <= 0) EXIT ! end of data?
6264  END DO
6265 
6266  ! release read buffer
6271 
6272  indab=0
6273  DO i=1,ntgb
6274  IF(globalparlabelindex(2,i) > 0) THEN
6275  IF(newcounter(i) >= mreqenf .OR. globalparcounts(i) >= iteren) THEN
6276  indab=indab+1
6277  globalparlabelindex(2,i)=indab ! variable, used in matrix (active)
6278  ELSE
6279  globalparlabelindex(2,i)=-3 ! fixed (iterated entries cut), not used in matrix (not active)
6280  END IF
6281  END IF
6282  END DO
6283  globalparheader(-6)=indab ! counted variable
6284  nvgb=indab ! nr of variable parameters
6285  WRITE(lunlog,*) 'LOOP1:',nvgb, ' is number NVGB of variable parameters'
6286  CALL mpdealloc(newcounter)
6287 
6288 END SUBROUTINE loop1i
6289 
6300 
6301 SUBROUTINE loop2
6302  USE mpmod
6303  USE mpdalc
6304 
6305  IMPLICIT NONE
6306  REAL(mps) :: chin2
6307  REAL(mps) :: chin3
6308  REAL(mps) :: cpr
6309  REAL(mps) :: fsum
6310  REAL(mps) :: gbc
6311  REAL(mps) :: gbu
6312  REAL(mpr8) :: glder
6313  INTEGER(mpi) :: i
6314  INTEGER(mpi) :: ia
6315  INTEGER(mpi) :: ib
6316  INTEGER(mpi) :: ibuf
6317  INTEGER(mpi) :: icblst
6318  INTEGER(mpi) :: icboff
6319  INTEGER(mpi) :: icgb
6320  INTEGER(mpi) :: iext
6321  INTEGER(mpi) :: ihis
6322  INTEGER(mpi) :: ij
6323  INTEGER(mpi) :: ij1
6324  INTEGER(mpi) :: ijn
6325  INTEGER(mpi) :: inder
6326  INTEGER(mpi) :: ioff
6327  INTEGER(mpi) :: ipoff
6328  INTEGER(mpi) :: iproc
6329  INTEGER(mpi) :: irecmm
6330  INTEGER(mpi) :: isfrst
6331  INTEGER(mpi) :: islast
6332  INTEGER(mpi) :: ist
6333  INTEGER(mpi) :: itgbi
6334  INTEGER(mpi) :: itgbij
6335  INTEGER(mpi) :: itgbik
6336  INTEGER(mpi) :: ivgbij
6337  INTEGER(mpi) :: ivgbik
6338  INTEGER(mpi) :: ivpgrp
6339  INTEGER(mpi) :: j
6340  INTEGER(mpi) :: ja
6341  INTEGER(mpi) :: jb
6342  INTEGER(mpi) :: jext
6343  INTEGER(mpi) :: jcgb
6344  INTEGER(mpi) :: jsp
6345  INTEGER(mpi) :: joff
6346  INTEGER(mpi) :: k
6347  INTEGER(mpi) :: kfile
6348  INTEGER(mpi) :: l
6349  INTEGER(mpi) :: label
6350  INTEGER(mpi) :: labelf
6351  INTEGER(mpi) :: labell
6352  INTEGER(mpi) :: lvpgrp
6353  INTEGER(mpi) :: lu
6354  INTEGER(mpi) :: lun
6355  INTEGER(mpi) :: maeqnf
6356  INTEGER(mpi) :: nall
6357  INTEGER(mpi) :: naeqna
6358  INTEGER(mpi) :: naeqnf
6359  INTEGER(mpi) :: naeqng
6360  INTEGER(mpi) :: nc31
6361  INTEGER(mpi) :: ncachd
6362  INTEGER(mpi) :: ncachi
6363  INTEGER(mpi) :: ncachr
6364  INTEGER(mpi) :: ncon
6365  INTEGER(mpi) :: nda
6366  INTEGER(mpi) :: ndf
6367  INTEGER(mpi) :: ndfmax
6368  INTEGER(mpi) :: nfixed
6369  INTEGER(mpi) :: nggd
6370  INTEGER(mpi) :: nggi
6371  INTEGER(mpi) :: nmatmo
6372  INTEGER(mpi) :: noff
6373  INTEGER(mpi) :: npar
6374  INTEGER(mpi) :: nparmx
6375  INTEGER(mpi) :: nr
6376  INTEGER(mpi) :: nrece
6377  INTEGER(mpi) :: nrecf
6378  INTEGER(mpi) :: nrecmm
6379  INTEGER(mpi) :: nst
6380  INTEGER(mpi) :: nwrd
6381  INTEGER(mpi) :: inone
6382  INTEGER(mpi) :: inc
6383  REAL(mps) :: wgh
6384  REAL(mps) :: wolfc3
6385  REAL(mps) :: wrec
6386  REAL(mps) :: chindl
6387 
6388  REAL(mpd)::dstat(3)
6389  REAL(mpd)::rerr
6390  INTEGER(mpl):: noff8
6391  INTEGER(mpl):: ndimbi
6392  INTEGER(mpl):: ndimsa(4)
6393  INTEGER(mpl):: ndgn
6394  INTEGER(mpl):: matsiz(2)
6395  INTEGER(mpl):: matwords
6396  INTEGER(mpl):: length
6397  INTEGER(mpl):: rows
6398  INTEGER(mpl):: cols
6399  INTEGER(mpl), PARAMETER :: two=2
6400  INTEGER(mpi) :: maxGlobalPar = 0
6401  INTEGER(mpi) :: maxLocalPar = 0
6402  INTEGER(mpi) :: maxEquations = 0
6403 
6404  INTERFACE ! needed for assumed-shape dummy arguments
6405  SUBROUTINE ndbits(npgrp,ndims,nsparr,ihst)
6406  USE mpdef
6407  INTEGER(mpi), DIMENSION(:), INTENT(IN) :: npgrp
6408  INTEGER(mpl), DIMENSION(4), INTENT(OUT) :: ndims
6409  INTEGER(mpl), DIMENSION(:,:), INTENT(OUT) :: nsparr
6410  INTEGER(mpi), INTENT(IN) :: ihst
6411  END SUBROUTINE ndbits
6412  SUBROUTINE ckbits(npgrp,ndims)
6413  USE mpdef
6414  INTEGER(mpi), DIMENSION(:), INTENT(IN) :: npgrp
6415  INTEGER(mpl), DIMENSION(4), INTENT(OUT) :: ndims
6416  END SUBROUTINE ckbits
6417  SUBROUTINE spbits(npgrp,nsparr,nsparc)
6418  USE mpdef
6419  INTEGER(mpi), DIMENSION(:), INTENT(IN) :: npgrp
6420  INTEGER(mpl), DIMENSION(:,:), INTENT(IN) :: nsparr
6421  INTEGER(mpi), DIMENSION(:), INTENT(OUT) :: nsparc
6422  END SUBROUTINE spbits
6423  SUBROUTINE gpbmap(npgrp,npair)
6424  USE mpdef
6425  INTEGER(mpi), DIMENSION(:,:), INTENT(IN) :: npgrp
6426  INTEGER(mpi), DIMENSION(:), INTENT(OUT) :: npair
6427  END SUBROUTINE gpbmap
6428  END INTERFACE
6429 
6430  SAVE
6431 
6432  !$ INTEGER(mpi) :: OMP_GET_THREAD_NUM
6433 
6434  isfrst(ibuf)=readbufferpointer(ibuf)+1
6435  islast(ibuf)=readbufferdatai(readbufferpointer(ibuf))
6436  inder(i)=readbufferdatai(i)
6437  glder(i)=readbufferdatad(i)
6438  ! ...
6439  WRITE(lunlog,*) ' '
6440  WRITE(lunlog,*) 'LOOP2: starting'
6441  CALL mstart('LOOP2')
6442 
6443  ! two subarrays to get the global parameter indices, used in an event
6444  length=nvgb
6445  CALL mpalloc(globalindexusage,length,'global index')
6446  CALL mpalloc(backindexusage,length,'back index')
6447  backindexusage=0
6448  CALL mpalloc(globalindexranges,length,'global index ranges')
6450 
6451  ! prepare constraints - determine number of constraints NCGB
6452  ! - sort and split into blocks
6453  ! - update globalIndexRanges
6454  CALL prpcon
6455 
6456  IF (metsol == 3.AND.icelim <= 0) THEN
6457  ! decomposition: enforce elimination
6458  icelim=1
6459  WRITE(lunlog,*) ' Elimination for constraints enforced for solution by decomposition!'
6460  END IF
6461  IF (icelim > 0) THEN ! elimination
6462  nagb=nvgb ! total number of parameters
6463  napgrp=nvpgrp ! total number of parameter groups
6464  nfgb=nvgb-ncgb ! number of fit parameters
6465  nprecond(1)=0 ! number of constraints for preconditioner
6466  nprecond(2)=nfgb ! matrix size for preconditioner
6467  ELSE ! Lagrange multipliers
6468  nagb=nvgb+ncgb ! total number of parameters
6469  napgrp=nvpgrp+ncgb ! total number of parameter groups
6470  nfgb=nagb ! number of fit parameters
6471  nprecond(1)=ncgb ! number of constraints for preconditioner
6472  nprecond(2)=nvgb ! matrix size for preconditioner
6473  ENDIF
6474  noff8=int(nagb,mpl)*int(nagb-1,mpl)/2
6475 
6476  ! all (variable) parameter groups
6477  length=napgrp+1
6478  CALL mpalloc(globalallindexgroups,length,'all parameter groups, 1. index')
6480  ivpgrp=0
6481  lvpgrp=-1
6482  DO i=1,ntgb
6483  ij=globalparlabelindex(2,i)
6484  IF (ij <= 0) cycle ! variable ?
6485  IF (globalparlabelindex(4,i) /= lvpgrp) THEN
6486  ivpgrp=ivpgrp+1
6487  globalallindexgroups(ivpgrp)=ij ! first index
6488  lvpgrp=globalparlabelindex(4,i)
6489  END IF
6490  END DO
6491  ! Lagrange multipliers
6492  IF (napgrp > nvpgrp) THEN
6493  DO jcgb=1, ncgb
6494  ivpgrp=ivpgrp+1
6495  globalallindexgroups(ivpgrp)=nvgb+jcgb
6496  END DO
6497  END IF
6499  ! from all (variable) parameters to group
6500  length=nagb
6501  CALL mpalloc(globalallpartogroup,length,'translation table all (var) par -> group')
6503  DO i=1,napgrp
6505  globalallpartogroup(j)=i
6506  END DO
6507  END DO
6508  IF (icheck > 2) THEN
6509  print *
6510  print *, ' Variable parameter groups ', nvpgrp
6511  DO i=1,nvpgrp
6514  globalparlabelindex(1,itgbi)
6515  END DO
6516  print *
6517  END IF
6518 
6519  ! read all data files and add all variable index pairs -------------
6520 
6521  IF (icheck > 1) CALL clbmap(ntpgrp)
6522 
6523  IF(matsto == 2) THEN
6524  CALL clbits(napgrp,mreqpe,mhispe,msngpe,mextnd,ndimbi,nspc) ! get dimension for bit storage, encoding, precision info
6525  END IF
6526 
6527  IF (imonit /= 0) THEN
6528  length=ntgb
6529  CALL mpalloc(measindex,length,'measurement counter/index')
6530  measindex=0
6531  CALL mpalloc(measres,length,'measurement resolution')
6532  measres=0.0_mps
6533  lunmon=9
6534  CALL mvopen(lunmon,'millepede.mon')
6535  ENDIF
6536 
6537  ! reading events===reading events===reading events===reading events=
6538  nrece =0 ! 'empty' records (no variable global parameters)
6539  nrecf =0 ! records with fixed global parameters
6540  naeqng=0 ! count number of equations (with global der.)
6541  naeqnf=0 ! count number of equations ( " , fixed)
6542  naeqna=0 ! all
6543  WRITE(lunlog,*) 'LOOP2: start event reading'
6544  ! monitoring for sparse matrix?
6545  irecmm=0
6546  IF (matsto == 2.AND.matmon /= 0) THEN
6547  nmatmo=0
6548  IF (matmon > 0) THEN
6549  nrecmm=matmon
6550  ELSE
6551  nrecmm=1
6552  END IF
6553  END IF
6554  DO k=1,3
6555  dstat(k)=0.0_mpd
6556  END DO
6557  ! define read buffer
6558  nc31=ncache/(31*mthrdr) ! split read cache 1 : 10 : 10*2 for pointers, ints, floats
6559  nwrd=nc31+1
6560  length=nwrd*mthrdr
6561  CALL mpalloc(readbufferpointer,length,'read buffer, pointer')
6562  nwrd=nc31*10+2+ndimbuf
6563  length=nwrd*mthrdr
6564  CALL mpalloc(readbufferdatai,length,'read buffer, integer')
6565  CALL mpalloc(readbufferdatad,length,'read buffer, real')
6566  ! to read (old) float binary files
6567  length=(ndimbuf+2)*