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