Millepede-II  V04-00-00_preview
 All Classes Files Functions Variables Enumerator Pages
pede.f90
Go to the documentation of this file.
1 
2 ! Code converted using TO_F90 by Alan Miller
3 ! Date: 2012-03-16 Time: 11:06:00
4 
87 
194 
348 
350 PROGRAM mptwo
351  USE mpmod
352  USE mpdalc
353  USE mptest1, ONLY: nplan,del,dvd
354  USE mptest2, ONLY: nlyr,nmx,nmy,sdevx,sdevy
355 
356  IMPLICIT NONE
357  REAL :: andf
358  REAL :: c2ndf
359  REAL :: deltat
360  REAL :: diff
361  REAL :: err
362  REAL :: gbu
363  REAL :: gmati
364  REAL :: rej
365  REAL :: rloop1
366  REAL :: rloop2
367  REAL :: rstext
368  REAL :: secnd
369  REAL :: rst
370  REAL ::rstp
371  REAL, DIMENSION(2) :: ta
372  INTEGER :: i
373  INTEGER :: ii
374  INTEGER :: ix
375  INTEGER :: ixv
376  INTEGER :: iy
377  INTEGER :: k
378  INTEGER :: kfl
379  INTEGER :: lun
380  INTEGER :: minut
381  INTEGER :: nhour
382  INTEGER :: nmxy
383  INTEGER :: nrc
384  INTEGER :: nsecnd
385  INTEGER :: ntot
386  INTEGER :: ntsec
387 
388  CHARACTER (LEN=24) :: chdate
389  CHARACTER (LEN=24) :: chost
390 
391  INTEGER(KIND=large) :: rows
392  INTEGER(KIND=large) :: cols
393 
394  DOUBLE PRECISION :: sums(9)
395  !$ INTEGER :: OMP_GET_NUM_PROCS,OMP_GET_MAX_THREADS
396  !$ INTEGER :: MXTHRD
397  !$ INTEGER :: NPROC
398 
399  SAVE
400  ! ...
401  CALL etime(ta,rstp)
402  CALL fdate(chdate)
403 
404  ! millepede.log file
405  lunlog=8
406  lvllog=1
407  CALL mvopen(lunlog,'millepede.log')
408  CALL getenv('HOSTNAME',chost)
409  IF (chost(1:1) == ' ') CALL getenv('HOST',chost)
410  WRITE(*,*) '($Rev: 92 $)'
411  !$ WRITE(*,*) 'using OpenMP (TM)'
412  !#ifdef __GFORTRAN__
413  ! WRITE(*,111) __GNUC__ , __GNUC_MINOR__ , __GNUC_PATCHLEVEL__
414  !111 FORMAT(' compiled with gcc ',i0,'.',i0,'.',i0)
415  !#endif
416  WRITE(*,*) ' '
417  WRITE(*,*) ' < Millepede II-P starting ... ',chdate
418  WRITE(*,*) ' ',chost
419  WRITE(*,*) ' '
420 
421  WRITE(8,*) ' '
422  WRITE(8,*) 'Log-file Millepede II-P ', chdate
423  WRITE(8,*) ' ', chost
424  ! CALL MEGDEB ! debug switched on
425 
426  ! read command line and text files
427 
428  CALL filetc ! command line and steering file analysis
429  CALL filetx ! read text files
430  lvllog=mprint ! export print level
431  IF (memdbg > 0) printflagalloc=1 ! debug memory management
432  !$ WRITE(*,*)
433  !$ NPROC=1
434  !$ MXTHRD=1
435  !$ NPROC=OMP_GET_NUM_PROCS() ! number of processors available
436  !$ CALL OMP_SET_NUM_THREADS(MTHRD) ! set max number of threads to MTHRD
437  !$ MXTHRD=OMP_GET_MAX_THREADS() ! get max number of threads back
438  !$ WRITE(*,*) 'Number of processors available: ', NPROC
439  !$ WRITE(*,*) 'Maximum number of OpenMP threads: ', MXTHRD
440  !$ WRITE(*,*) 'Number of threads for processing: ', MTHRD
441  !$ IF (MXREC.GT.0) MTHRDR=1 ! to get allways the same MXREC records
442  !$ WRITE(*,*) 'Number of threads for reading: ', MTHRDR
443  !$POMP INST INIT ! start profiling with ompP
444  IF (ncache < 0) THEN
445  ncache=25000000*mthrd ! default cache size (100 MB per thread)
446  ENDIF
447  rows=6; cols=mthrdr
448  CALL mpalloc(readbufferinfo,rows,cols,'read buffer header')
449  ! histogram file
450  lun=7
451  CALL mvopen(lun,'millepede.his')
452  CALL hmplun(lun) ! unit for histograms
453  CALL gmplun(lun) ! unit for xy data
454 
455  ! debugging
456  IF(nrecpr /= 0.OR.nrecp2 /= 0) THEN
457  CALL mvopen(1,'mpdebug.txt')
458  END IF
459 
460  CALL etime(ta,rstext)
461  times(0)=rstext-rstp ! time for text processing
462 
463  ! preparation of data sub-arrays
464 
465  CALL loop1
466  CALL etime(ta,rloop1)
467  times(1)=rloop1-rstext ! time for LOOP1
468 
469  CALL loop2
470  IF(chicut /= 0.0) THEN
471  WRITE(8,*) 'Chi square cut equiv 3 st.dev applied ...'
472  WRITE(8,*) ' in first iteration with factor',chicut
473  WRITE(8,*) ' in second iteration with factor',chirem
474  WRITE(8,*) ' (reduced by sqrt in next iterations)'
475  END IF
476  IF(lhuber /= 0) THEN
477  WRITE(8,*) 'Down-weighting of outliers in', lhuber,' iterations'
478  WRITE(8,*) 'Cut on downweight fraction',dwcut
479  END IF
480 
481  CALL etime(ta,rloop2)
482  times(2)=rloop2-rloop1 ! time for LOOP2
483 
484  ! use different solution methods
485 
486  CALL mstart('Iteration') ! Solution module starting
487 
488  ! CALL XEXECS ! all methods
489 
490  CALL xloopn ! all methods
491 
492  ! IF(NREJEC(1)+NREJEC(2)+NREJEC(3).NE.0) THEN
493  ! WRITE(*,*) 'Data rejected in last loop: ',NREJEC(1),
494  ! + ' (Ndf=0) ',NREJEC(2),' (huge) ',NREJEC(3),' (large)'
495  ! END IF
496 
497 
498 
499  ! ------------------------------------------------------------------
500 
501  IF(nloopn > 2.AND.nhistp /= 0) THEN ! last iteration
502  CALL hmprnt(3) ! scaled residual of single measurement (with global deriv.)
503  CALL hmprnt(12) ! scaled residual of single measurement (no global deriv.)
504  CALL hmprnt(4) ! chi^2/Ndf
505  END IF
506  IF(nloopn > 2) THEN
507  CALL hmpwrt(3)
508  CALL hmpwrt(12)
509  CALL hmpwrt(4)
510  CALL gmpwrt(4) ! location, dispersion (res.) as a function of record nr
511  IF (nloopn <= lfitnp) THEN
512  CALL hmpwrt(13)
513  CALL hmpwrt(14)
514  CALL gmpwrt(5)
515  END IF
516  END IF
517  IF(nhistp /= 0) THEN
518  CALL gmprnt(1)
519  CALL gmprnt(2)
520  END IF
521  CALL gmpwrt(1) ! output of xy data
522  CALL gmpwrt(2) ! output of xy data
523  ! 'track quality' per binary file
524  IF (nfilb > 1) THEN
525  CALL gmpdef(6,1,'log10(#records) vs file number')
526  CALL gmpdef(7,1,'final rejection fraction vs file number')
527  CALL gmpdef(8,1, &
528  'final <Chi^2/Ndf> from accepted local fits vs file number')
529  CALL gmpdef(9,1, '<Ndf> from accepted local fits vs file number')
530 
531  DO i=1,nfilb
532  kfl=kfd(2,i)
533  nrc=-kfd(1,i)
534  IF (nrc > 0) THEN
535  rej=float(nrc-jfd(kfl))/float(nrc)
536  CALL gmpxy(6,float(kfl),alog10(float(nrc))) ! log10(#records) vs file
537  CALL gmpxy(7,float(kfl),rej) ! rejection fraction vs file
538  END IF
539  IF (jfd(kfl) > 0) THEN
540  c2ndf=cfd(kfl)/float(jfd(kfl))
541  CALL gmpxy(8,float(kfl),c2ndf) ! <Chi2/NDF> vs file
542  andf=float(dfd(kfl))/float(jfd(kfl))
543  CALL gmpxy(9,float(kfl),andf) ! <NDF> vs file
544  END IF
545  END DO
546  IF(nhistp /= 0) THEN
547  CALL gmprnt(6)
548  CALL gmprnt(7)
549  CALL gmprnt(8)
550  CALL gmprnt(9)
551  END IF
552  CALL gmpwrt(6) ! output of xy data
553  CALL gmpwrt(7) ! output of xy data
554  CALL gmpwrt(8) ! output of xy data
555  CALL gmpwrt(9) ! output of xy data
556  END IF
557 
558  IF(ictest == 1) THEN
559  WRITE(*,*) ' '
560  WRITE(*,*) 'Misalignment test wire chamber'
561  WRITE(*,*) ' '
562 
563  CALL hmpdef( 9,-0.0015,+0.0015,'True - fitted displacement')
564  CALL hmpdef(10,-0.0015,+0.0015,'True - fitted Vdrift')
565  DO i=1,4
566  sums(i)=0.0d0
567  END DO
568  DO i=1,nplan
569  diff=sngl(-del(i)-globalparameter(i))
570  sums(1)=sums(1)+diff
571  sums(2)=sums(2)+diff*diff
572  diff=sngl(-dvd(i)-globalparameter(100+i))
573  sums(3)=sums(3)+diff
574  sums(4)=sums(4)+diff*diff
575  END DO
576  sums(1)=0.01d0*sums(1)
577  sums(2)=sqrt(0.01d0*sums(2))
578  sums(3)=0.01d0*sums(3)
579  sums(4)=sqrt(0.01d0*sums(4))
580  WRITE(*,143) 'Parameters 1 - 100: mean =',sums(1), 'rms =',sums(2)
581  WRITE(*,143) 'Parameters 101 - 200: mean =',sums(3), 'rms =',sums(4)
582 143 format(6x,a28,f9.6,3x,a5,f9.6)
583  WRITE(*,*) ' '
584  WRITE(*,*) ' '
585  WRITE(*,*) ' I '
586  WRITE(*,*) ' --- '
587  DO i=1,100
588  WRITE(*,102) i,-del(i),globalparameter(i),-del(i)-globalparameter(i), &
589  -dvd(i),globalparameter(100+i),-dvd(i)-globalparameter(100+i)
590  diff=sngl(-del(i)-globalparameter(i))
591  CALL hmpent( 9,diff)
592  diff=sngl(-dvd(i)-globalparameter(100+i))
593  CALL hmpent(10,diff)
594  END DO
595  IF(nhistp /= 0) THEN
596  CALL hmprnt( 9)
597  CALL hmprnt(10)
598  END IF
599  CALL hmpwrt( 9)
600  CALL hmpwrt(10)
601  END IF
602  IF(ictest > 1) THEN
603  WRITE(*,*) ' '
604  WRITE(*,*) 'Misalignment test Si tracker'
605  WRITE(*,*) ' '
606 
607  CALL hmpdef( 9,-0.0025,+0.0025,'True - fitted displacement X')
608  CALL hmpdef(10,-0.025,+0.025,'True - fitted displacement Y')
609  DO i=1,9
610  sums(i)=0.0d0
611  END DO
612  nmxy=nmx*nmy
613  ix=0
614  iy=ntot
615  DO i=1,nlyr
616  DO k=1,nmxy
617  ix=ix+1
618  diff=sngl(-sdevx((i-1)*nmxy+k)-globalparameter(ix))
619  sums(1)=sums(1)+1.0d0
620  sums(2)=sums(2)+diff
621  sums(3)=sums(3)+diff*diff
622  ixv=globalparlabelindex(2,ix)
623  IF (ixv > 0.AND.metsol == 1.OR.metsol == 2) THEN
624  ii=(ixv*ixv+ixv)/2
625  gmati=sngl(globalmatd(ii))
626  err=sqrt(abs(gmati))
627  diff=diff/err
628  sums(7)=sums(7)+1.0d0
629  sums(8)=sums(8)+diff
630  sums(9)=sums(9)+diff*diff
631  END IF
632  END DO
633  IF (mod(i,3) == 1) THEN
634  DO k=1,nmxy
635  iy=iy+1
636  diff=-sngl(sdevy((i-1)*nmxy+k)-globalparameter(iy))
637  sums(4)=sums(4)+1.0d0
638  sums(5)=sums(5)+diff
639  sums(6)=sums(6)+diff*diff
640  ixv=globalparlabelindex(2,iy)
641  IF (ixv > 0.AND.metsol == 1.OR.metsol == 2) THEN
642  ii=(ixv*ixv+ixv)/2
643  gmati=sngl(globalmatd(ii))
644  err=sqrt(abs(gmati))
645  diff=diff/err
646  sums(7)=sums(7)+1.0d0
647  sums(8)=sums(8)+diff
648  sums(9)=sums(9)+diff*diff
649  END IF
650  END DO
651  END IF
652  END DO
653  sums(2)=sums(2)/sums(1)
654  sums(3)=sqrt(sums(3)/sums(1))
655  sums(5)=sums(5)/sums(4)
656  sums(6)=sqrt(sums(6)/sums(4))
657  WRITE(*,143) 'Parameters 1 - 500: mean =',sums(2), 'rms =',sums(3)
658  WRITE(*,143) 'Parameters 501 - 700: mean =',sums(5), 'rms =',sums(6)
659  IF (sums(7) > 0.5d0) THEN
660  sums(8)=sums(8)/sums(7)
661  sums(9)=sqrt(sums(9)/sums(7))
662  WRITE(*,143) 'Parameter pulls, all: mean =',sums(8), 'rms =',sums(9)
663  END IF
664  WRITE(*,*) ' '
665  WRITE(*,*) ' '
666  WRITE(*,*) ' I '
667  WRITE(*,*) ' --- '
668  ix=0
669  iy=ntot
670  DO i=1,nlyr
671  DO k=1,nmxy
672  ix=ix+1
673  diff=sngl(-sdevx((i-1)*nmxy+k)-globalparameter(ix))
674  CALL hmpent( 9,diff)
675  WRITE(*,102) ix,-sdevx((i-1)*nmxy+k),globalparameter(ix),-diff
676  END DO
677  END DO
678  DO i=1,nlyr
679  IF (mod(i,3) == 1) THEN
680  DO k=1,nmxy
681  iy=iy+1
682  diff=sngl(-sdevy((i-1)*nmxy+k)-globalparameter(iy))
683  CALL hmpent(10,diff)
684  WRITE(*,102) iy,-sdevy((i-1)*nmxy+k),globalparameter(iy),-diff
685  END DO
686  END IF
687  END DO
688  IF(nhistp /= 0) THEN
689  CALL hmprnt( 9)
690  CALL hmprnt(10)
691  END IF
692  CALL hmpwrt( 9)
693  CALL hmpwrt(10)
694  END IF
695 
696  IF(nrec1+nrec2 > 0) THEN
697  WRITE(8,*) ' '
698  IF(nrec1 > 0) THEN
699  WRITE(8,*) 'Record',nrec1,' has largest residual:',value1
700  END IF
701  IF(nrec2 > 0) THEN
702  WRITE(8,*) 'Record',nrec2,' has largest Chi^2/Ndf:',value2
703  END IF
704  END IF
705  IF(nrec3 < maxi4) THEN
706  WRITE(8,*) 'Record',nrec3, ' is first with error (rank deficit/NaN)'
707  END IF
708  WRITE(8,*) ' '
709  WRITE(8,*) 'In total 2 +',nloopn,' loops through the data files'
710  IF (mnrsit > 0) THEN
711  WRITE(8,*) ' '
712  WRITE(8,*) 'In total ',mnrsit,' internal MINRES iterations'
713  END IF
714 
715  WRITE(8,103) times(0),times(1),times(2),times(4),times(7), &
716  times(5),times(8),times(3),times(6)
717  CALL etime(ta,rst)
718  deltat=rst-rstp
719  ntsec=nint(deltat)
720  CALL sechms(deltat,nhour,minut,secnd)
721  nsecnd=nint(secnd) ! round
722  WRITE(8,*) 'Total time =',ntsec,' seconds =',nhour,' h',minut, &
723  ' m',nsecnd,' seconds'
724  CALL fdate(chdate)
725  WRITE(8,*) 'end ', chdate
726  gbu=4.02e-9*float(maxwordsalloc) ! GB used
727  WRITE(8,*) ' '
728  WRITE(8,105) gbu
729 
730  ! Rejects ----------------------------------------------------------
731 
732  IF(nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) /= 0) THEN
733  WRITE(8,*) ' '
734  WRITE(8,*) 'Data rejected in last iteration: '
735  WRITE(8,*) ' ', &
736  nrejec(0), ' (rank deficit/NaN) ',nrejec(1),' (Ndf=0) ', &
737  nrejec(2), ' (huge) ',nrejec(3),' (large)'
738  WRITE(8,*) ' '
739  END IF
740  CALL explfc(8)
741 
742  WRITE(*,*) ' '
743  WRITE(*,*) ' < Millepede II-P ending ... ', chdate ! with exit code',ITEXIT,' >'
744  WRITE(*,*) ' '
745  gbu=4.0e-9*float(maxwordsalloc) ! GB used
746  WRITE(*,105) gbu
747  WRITE(*,*) ' '
748 
749 102 format(2x,i4,2x,3f10.5,2x,3f10.5)
750 103 format(' Times [in sec] for text processing',f12.3/ &
751  ' LOOP1',f12.3/ &
752  ' LOOP2',f12.3/ &
753  ' func. value ',f12.3,' *',f4.0/ &
754  ' func. value, global matrix, solution',f12.3,' *',f4.0/ &
755  ' new solution',f12.3,' *',f4.0/)
756 105 format(' Peak dynamic memory allocation: ',f11.6,' GB')
757 END PROGRAM mptwo ! Mille
758 
759 
766 
767 SUBROUTINE solglo(ivgbi)
768  USE mpmod
769 
770  IMPLICIT NONE
771  REAL :: dpa
772  REAL :: err
773  REAL :: gcor2
774  REAL :: par
775  INTEGER :: iph
776  INTEGER :: istop
777  INTEGER :: itgbi
778  INTEGER :: itgbl
779  INTEGER :: itn
780  INTEGER :: itnlim
781  INTEGER :: nout
782 
783  INTEGER, INTENT(IN) :: ivgbi
784 
785  DOUBLE PRECISION :: shift
786  DOUBLE PRECISION :: rtol
787  DOUBLE PRECISION :: anorm
788  DOUBLE PRECISION :: acond
789  DOUBLE PRECISION :: rnorm
790  DOUBLE PRECISION :: ynorm
791  DOUBLE PRECISION :: gmati
792  DOUBLE PRECISION :: diag
793  INTEGER(KIND=large) :: ijadd
794  INTEGER(KIND=large) :: jk
795  INTEGER(KIND=large) :: ii
796  LOGICAL :: checka
797  EXTERNAL avprod, mcsolv, mvsolv
798  SAVE
799  DATA iph/0/
800  ! ...
801  IF(iph == 0) THEN
802  iph=1
803  WRITE(*,101)
804  END IF
805  itgbi=globalparvartototal(ivgbi)
806  itgbl=globalparlabelindex(1,itgbi)
807 
808  globalvector=0.0d0 ! reset rhs vector IGVEC
809  globalvector(ivgbi)=1.0d0
810 
811  ! NOUT =6
812  nout =0
813  itnlim=200
814  shift =0.0d0
815  rtol = mrestl ! from steering
816  checka=.false.
817 
818  IF(mbandw == 0) THEN ! default preconditioner
819  CALL minres(nagb,globalvector, &
820  workspaced(1), workspaced(1+nagb), workspaced(1+2*nagb), &
821  workspaced(1+3*nagb),workspaced(1+4*nagb),workspaced(1+5*nagb), &
822  globalcorrections,workspaceminres, avprod, mcsolv, checka ,.true. , shift, &
823  nout , itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm)
824  ELSE IF(mbandw > 0) THEN ! band matrix preconditioner
825  CALL minres(nagb,globalvector, &
826  workspaced(1), workspaced(1+nagb), workspaced(1+2*nagb), &
827  workspaced(1+3*nagb),workspaced(1+4*nagb),workspaced(1+5*nagb), &
828  globalcorrections,workspaceminres, avprod, mvsolv, checka ,.true. , shift, &
829  nout, itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm)
830  ELSE
831  CALL minres(nagb,globalvector, &
832  workspaced(1), workspaced(1+nagb), workspaced(1+2*nagb), &
833  workspaced(1+3*nagb),workspaced(1+4*nagb),workspaced(1+5*nagb), &
834  globalcorrections,workspaceminres, avprod, mvsolv, checka ,.false., shift, &
835  nout, itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm)
836  END IF
837 
838  ! subroutine MINRES( n, b, r1, r2, v, w, w1, w2, x, y,
839  ! $ AVPROD, Msolve, .TRUE. , .FALSE. , SHIFT,
840  ! $ NOUT , ITNLIM, rtol,
841  ! $ ISTOP, ITN, Anorm, Acond, rnorm, ynorm )
842 
843  par=sngl(globalparameter(itgbi))
844  dpa=par-globalparstart(itgbi)
845  gmati=globalcorrections(ivgbi)
846  err=sqrt(abs(sngl(gmati)))
847  IF(gmati < 0.0d0) err=-err
848  IF(matsto == 1) THEN ! normal matrix ! ???
849  ii=ivgbi
850  jk=(ii*ii+ii)/2
851  ELSE IF(matsto == 2) THEN ! sparse matrix
852  jk=ijadd(ivgbi,ivgbi)
853  END IF
854  IF (jk > 0) THEN
855  diag=globalmatd(jk)
856  ELSE
857  diag=dble(globalmatf(-jk))
858  END IF
859  gcor2=sngl(1.0d0-1.0d0/(gmati*diag)) ! global correlation (squared)
860  WRITE(*,102) itgbl,par,globalparpresigma(itgbi),dpa,err,gcor2,itn
861 101 format(1x,' label parameter presigma differ', &
862  ' Error gcor^2 iit'/ 1x,'---------',2x,5('-----------'),2x,'----')
863 102 format(i10,2x,4g12.4,f7.4,i6,i4)
864 END SUBROUTINE solglo
865 
867 SUBROUTINE addcst
868  USE mpmod
869 
870  IMPLICIT NONE
871  REAL :: climit
872  REAL :: factr
873  REAL :: hugeval
874  REAL :: sgm
875 
876  INTEGER :: i
877  INTEGER :: icgb
878  INTEGER :: irhs
879  INTEGER :: itgbi
880  INTEGER :: ivgb
881  INTEGER :: l
882  INTEGER :: label
883  INTEGER :: nop
884  INTEGER :: inone
885 
886  DOUBLE PRECISION :: rhs
887  DOUBLE PRECISION :: drhs(4)
888  INTEGER :: idrh (4)
889  SAVE
890  ! ...
891  nop=0
892  hugeval=1.0 ! E6
893  IF(lenconstraints == 0) return ! no constraints
894  climit=1.0e-5 ! limit for printout
895  irhs=0 ! number of values in DRHS(.), to be printed
896 
897  i=1
898  DO icgb=1,ncgb
899  rhs=listconstraints(i )%value ! right hand side
900  sgm=listconstraints(i+1)%value ! sigma parameter
901  i=i+2
902  DO
903  label=listconstraints(i)%label
904  factr=listconstraints(i)%value
905  itgbi=inone(label) ! -> ITGBI= index of parameter label
906  ivgb =globalparlabelindex(2,itgbi) ! -> variable-parameter index
907 
908  IF(icalcm == 1.AND.ivgb > 0) THEN
909  CALL mupdat(nvgb+icgb,ivgb,dble(hugeval*factr)) ! add to matrix
910  END IF
911 
912  rhs=rhs-factr*globalparameter(itgbi) ! reduce residuum
913  i=i+1
914  IF(i > lenconstraints) exit
915  IF(listconstraints(i)%label == 0) exit
916  END DO
917  IF(abs(rhs) > climit) THEN
918  irhs=irhs+1
919  idrh(irhs)=icgb
920  drhs(irhs)=rhs
921  nop=1
922  IF(irhs == 4) THEN
923  WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
924  irhs=0
925  END IF
926  END IF
927  IF(rhs == 0.0) rhs=1.0e-9
928  globalvector(nvgb+icgb)=rhs*hugeval
929  END DO
930 
931  IF(irhs /= 0) THEN
932  WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
933  END IF
934  IF(nop == 0) return
935  WRITE(*,102) ' Constraints: only equation values >', climit,' are printed'
936 101 format(' ',4(i4,g11.3))
937 102 format(a,g11.2,a)
938 END SUBROUTINE addcst
939 
943 
944 SUBROUTINE feasma
945  USE mpmod
946  USE mpdalc
947 
948  IMPLICIT NONE
949  REAL :: factr
950  REAL :: sgm
951  INTEGER :: i
952  INTEGER :: icgb
953  INTEGER :: ij
954  INTEGER :: itgbi
955  INTEGER :: ivgb
956  INTEGER :: j
957  INTEGER :: l
958  INTEGER :: label
959  INTEGER :: nrank
960  INTEGER :: inone
961 
962  DOUBLE PRECISION:: rhs
963  INTEGER(KIND=large):: length
964  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: vecconstraints
965  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: matconstraintst
966  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: auxvectord
967  INTEGER, DIMENSION(:), ALLOCATABLE :: auxvectori
968  SAVE
969  ! ...
970  IF(lenconstraints == 0) return ! no constraints
971 
972  length=(ncgb*ncgb+ncgb)/2
973  CALL mpalloc(matconsproduct, length, 'product matrix of constraints')
974  matconsproduct=0.0d0
975  length=ncgb
976  CALL mpalloc(vecconsresiduals, length, 'residuals of constraints')
977  CALL mpalloc(vecconstraints,length,'rhs vector') ! double precision vector
978  CALL mpalloc(auxvectori,length,'auxiliary array (I)') ! int aux 1
979  CALL mpalloc(auxvectord,length,'auxiliary array (D)') ! double aux 1
980  ! constraint matrix A and product A A^T (A is stored as transposed)
981  length = ncgb*nvgb
982  CALL mpalloc(matconstraintst,length,'transposed matrix of constraints')
983  matconstraintst=0.0d0
984 
985  i=1
986  DO icgb=1,ncgb
987  rhs=listconstraints(i )%value ! right hand side
988  sgm=listconstraints(i+1)%value ! sigma parameter
989  i=i+2
990  DO
991  label=listconstraints(i)%label
992  factr=listconstraints(i)%value
993  itgbi=inone(label) ! -> ITGBI= index of parameter label
994  ivgb =globalparlabelindex(2,itgbi) ! -> variable-parameter index
995  IF(ivgb > 0) matconstraintst(icgb+(ivgb-1)*ncgb)=factr ! matrix element
996  rhs=rhs-factr*globalparameter(itgbi) ! reduce residuum
997  i=i+1
998  IF(i > lenconstraints) exit
999  IF(listconstraints(i)%label == 0) exit
1000  END DO
1001  vecconstraints(icgb)=rhs ! constraint discrepancy
1002  END DO
1003 
1004  DO l=1,nvgb
1005  ij=0
1006  DO i=1,ncgb
1007  DO j=1,i
1008  ij=ij+1
1009  matconsproduct(ij)=matconsproduct(ij)+matconstraintst((l-1)*ncgb+i)*matconstraintst((l-1)*ncgb+j)
1010  END DO
1011  END DO
1012  END DO
1013 
1014  ! inversion of product matrix of constraints
1015  CALL sqminv(matconsproduct,vecconstraints,ncgb,nrank, auxvectord, auxvectori)
1016 
1017  nmiss1=ncgb-nrank
1018 
1019  WRITE(*,*) ' '
1020  WRITE(*,*) 'Rank of product matrix of constraints is',nrank, &
1021  ' for',ncgb,' constraint equations'
1022  WRITE(8,*) 'Rank of product matrix of constraints is',nrank, &
1023  ' for',ncgb,' constraint equations'
1024  IF(nrank < ncgb) THEN
1025  WRITE(*,*) 'Warning: insufficient constraint equations!'
1026  WRITE(8,*) 'Warning: insufficient constraint equations!'
1027  IF (iforce == 0) THEN
1028  isubit=1
1029  WRITE(*,*) ' --> enforcing SUBITO mode'
1030  WRITE(8,*) ' --> enforcing SUBITO mode'
1031  END IF
1032  END IF
1033 
1034  CALL mpdealloc(matconstraintst)
1035  CALL mpdealloc(auxvectord)
1036  CALL mpdealloc(auxvectori)
1037  CALL mpdealloc(vecconstraints)
1038 
1039  return
1040 END SUBROUTINE feasma ! matrix for feasible solution
1041 
1049 SUBROUTINE feasib(concut,iact)
1050  USE mpmod
1051  USE mpdalc
1052 
1053  IMPLICIT NONE
1054  REAL :: factr
1055  REAL :: sgm
1056  INTEGER :: i
1057  INTEGER :: icgb
1058  INTEGER :: iter
1059  INTEGER :: itgbi
1060  INTEGER :: ivgb
1061  INTEGER :: label
1062  INTEGER :: inone
1063 
1064  REAL, INTENT(IN) :: concut
1065  INTEGER, INTENT(OUT) :: iact
1066 
1067  DOUBLE PRECISION :: rhs
1068  DOUBLE PRECISION ::sum1
1069  DOUBLE PRECISION ::sum2
1070  DOUBLE PRECISION ::sum3
1071  INTEGER(KIND=large) :: length
1072 
1073  ! ! acts on subarray "0"
1074  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: auxvectord
1075  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: veccorrections
1076  SAVE
1077 
1078  iact=0
1079  IF(lenconstraints == 0) return ! no constraints
1080 
1081  DO iter=1,2
1082  vecconsresiduals=0.0d0
1083 
1084  ! calculate right constraint equation discrepancies
1085  i=1
1086  DO icgb=1,ncgb
1087  rhs=listconstraints(i )%value ! right hand side
1088  sgm=listconstraints(i+1)%value ! sigma parameter
1089  i=i+2
1090  DO
1091  label=listconstraints(i)%label
1092  factr=listconstraints(i)%value
1093  itgbi=inone(label) ! -> ITGBI= index of parameter label
1094  rhs=rhs-factr*globalparameter(itgbi) ! reduce residuum
1095  i=i+1
1096  IF(i > lenconstraints) exit
1097  IF(listconstraints(i)%label == 0) exit
1098  ENDDO
1099  vecconsresiduals(icgb)=rhs ! constraint discrepancy
1100  END DO
1101 
1102  ! constraint equation discrepancies -------------------------------
1103 
1104  sum1=0.0d0
1105  sum2=0.0d0
1106  sum3=0.0d0
1107  DO icgb=1,ncgb
1108  sum1=sum1+vecconsresiduals(icgb)**2
1109  sum2=sum2+abs(vecconsresiduals(icgb))
1110  sum3=max(sum3,abs(vecconsresiduals(icgb)))
1111  END DO
1112  sum1=sqrt(sum1/dfloat(ncgb))
1113  sum2=sum2/dfloat(ncgb)
1114 
1115  IF(iter == 1.AND.sum1 < concut) return ! do nothing if correction small
1116 
1117  IF(iter == 1.AND.ncgb <= 12) THEN
1118  WRITE(*,*) ' '
1119  WRITE(*,*) 'Constraint equation discrepancies:'
1120  WRITE(*,101) (icgb,vecconsresiduals(icgb),icgb=1,ncgb)
1121 101 format(4x,4(i5,g12.4))
1122  WRITE(*,103) concut
1123 103 format(10x,' Cut on rms value is',g8.1)
1124  END IF
1125 
1126  IF(iact == 0) THEN
1127  WRITE(*,*) ' '
1128  WRITE(*,*) 'Improve constraints'
1129  END IF
1130  iact=1
1131 
1132  WRITE(*,102) iter,sum1,sum2,sum3
1133 102 format(i6,' rms',g12.4,' avrg_abs',g12.4,' max_abs',g12.4)
1134 
1135  length=ncgb
1136  CALL mpalloc(auxvectord,length,'auxiliary vector (D)')
1137  length=nvgb
1138  CALL mpalloc(veccorrections,length,'constraint corrections')
1139  veccorrections=0.0d0
1140 
1141  ! multiply inverse matrix and constraint vector
1142  CALL dbsvx(matconsproduct,vecconsresiduals,auxvectord,ncgb)
1143 
1144  i=1
1145  DO icgb=1,ncgb
1146  rhs=listconstraints(i )%value ! right hand side
1147  sgm=listconstraints(i+1)%value ! sigma parameter
1148  i=i+2
1149  DO
1150  label=listconstraints(i)%label
1151  factr=listconstraints(i)%value
1152  itgbi=inone(label) ! -> ITGBI= index of parameter label
1153  ivgb =globalparlabelindex(2,itgbi) ! -> variable-parameter index
1154  IF(ivgb > 0) THEN
1155  veccorrections(ivgb)=veccorrections(ivgb)+auxvectord(icgb)*factr
1156  END IF
1157  i=i+1
1158  IF(i > lenconstraints) exit
1159  IF(listconstraints(i)%label == 0) exit
1160  ENDDO
1161  END DO
1162 
1163  DO i=1,nvgb ! add corrections
1164  itgbi=globalparvartototal(i)
1165  globalparameter(itgbi)=globalparameter(itgbi)+veccorrections(i)
1166  END DO
1167 
1168  CALL mpdealloc(veccorrections)
1169  CALL mpdealloc(auxvectord)
1170 
1171  END DO ! iteration 1 and 2
1172 
1173 END SUBROUTINE feasib ! make parameters feasible
1174 
1207 SUBROUTINE peread(more)
1208  USE mpmod
1209 
1210  IMPLICIT NONE
1211  INTEGER :: i
1212  INTEGER :: iact
1213  INTEGER :: ierrc
1214  INTEGER :: ierrf
1215  INTEGER :: inder
1216  INTEGER :: ioffp
1217  INTEGER :: ithr
1218  INTEGER :: jfile
1219  INTEGER :: jrec
1220  INTEGER :: k
1221  INTEGER :: kfile
1222  INTEGER :: l
1223  INTEGER :: lun
1224  INTEGER :: mpri
1225  INTEGER :: n
1226  INTEGER :: nact
1227  INTEGER :: nbuf
1228  INTEGER :: ndata
1229  INTEGER :: noff
1230  INTEGER :: npointer
1231  INTEGER :: npri
1232  INTEGER :: nr
1233  INTEGER :: nrc
1234  INTEGER :: nrd
1235  INTEGER :: nrpr
1236  INTEGER :: nthr
1237  INTEGER :: ntot
1238  INTEGER :: maxrecordsize
1239  INTEGER :: maxrecordfile
1240 
1241  INTEGER, INTENT(OUT) :: more
1242 
1243  LOGICAL :: lprint
1244  LOGICAL :: floop
1245  LOGICAL :: eof
1246  DOUBLE PRECISION :: ds0
1247  DOUBLE PRECISION :: ds1
1248  DOUBLE PRECISION :: ds2
1249  DOUBLE PRECISION :: dw
1250  !$ INTEGER :: OMP_GET_THREAD_NUM
1251  SAVE
1252 
1253  inder(i)=readbufferdatai(i)
1254 
1255  DATA lprint/.true./
1256  DATA floop/.true./
1257  DATA npri / 0 /, mpri / 1000 /
1258  ! ...
1259  IF(ifile == 0) THEN ! start/restart
1260  nrec=0
1261  ntot=0
1262  sumrecords=0
1263  skippedrecords=0
1264  numblocks=0
1265  minrecordsinblock=size(readbufferdatai)
1266  maxrecordsinblock=0
1267  readbufferinfo=0 ! reset management info
1268  nrpr=1
1269  nthr=mthrdr
1270  nact=0 ! active threads (have something still to read)
1271  DO k=1,nthr
1272  IF (ifile < nfilb) THEN
1273  ifile=ifile+1
1274  readbufferinfo(1,k)=ifile
1275  readbufferinfo(2,k)=nact
1276  nact=nact+1
1277  END IF
1278  END DO
1279  END IF
1280  npointer=size(readbufferpointer)/nact
1281  ndata=size(readbufferdatai)/nact
1282  more=-1
1283  DO k=1,nthr
1284  iact=readbufferinfo(2,k)
1285  readbufferinfo(4,k)=0 ! reset counter
1286  readbufferinfo(5,k)=iact*ndata ! reset offset
1287  END DO
1288  numblocks=numblocks+1 ! new block
1289 
1290  !$OMP PARALLEL &
1291  !$OMP DEFAULT(PRIVATE) &
1292  !$OMP SHARED(readBufferInfo,readBufferPointer,readBufferDataI,readBufferDataF, &
1293  !$OMP nPointer,nData,skippedRecords,ndimbuf,NTHR,NFILF,FLOOP, &
1294  !$OMP IFD,KFD,IFILE,NFILB,WFD,XFD) &
1295  !$OMP NUM_THREADS(NTHR)
1296 
1297  ithr=1
1298  !$ ITHR=OMP_GET_THREAD_NUM()+1 ! thread number
1299  jfile=readbufferinfo(1,ithr) ! file index
1300  iact =readbufferinfo(2,ithr) ! active thread number
1301  jrec =readbufferinfo(3,ithr) ! records read
1302  ioffp=iact*npointer
1303 
1304  files: DO WHILE (jfile > 0)
1305  kfile=kfd(2,jfile)
1306  records: DO
1307  nbuf=readbufferinfo(4,ithr)+1
1308  noff=readbufferinfo(5,ithr)+2 ! 2 header words per record
1309  nr=ndimbuf
1310  IF(kfile <= nfilf) THEN ! Fortran file
1311  lun=kfile+10
1312  READ(lun,iostat=ierrf) n,(readbufferdataf(noff+i),i=1,min(n/2,nr)),&
1313  (readbufferdatai(noff+i),i=1,min(n/2,nr))
1314  nr=n/2
1315  IF (ierrf < 0) rewind lun ! end-of-file
1316  eof=(ierrf /= 0)
1317  ELSE ! C file
1318  lun=kfile-nfilf
1319 #ifdef READ_C_FILES
1320  CALL readc(readbufferdataf(noff+1),readbufferdatai(noff+1),nr,lun,ierrc)
1321  n=nr+nr
1322 #else
1323  ierrc=0
1324 #endif
1325  eof=(ierrc <= 0.AND.ierrc /= -4) ! allow buffer overruns -> skip record
1326  END IF
1327  IF(eof) exit records ! end-of-files or error
1328 
1329  jrec=jrec+1
1330  readbufferinfo(3,ithr)=jrec
1331  IF(floop) THEN
1332  xfd(jfile)=max(xfd(jfile),n)
1333  IF(ithr == 1) THEN
1334  CALL hmplnt(1,n)
1335  IF(inder(noff+1) /= 0) CALL hmpent(8,float(inder(noff+1)))
1336  END IF
1337  END IF
1338 
1339  IF (nr <= ndimbuf) THEN
1340  readbufferinfo(4,ithr)=nbuf
1341  readbufferinfo(5,ithr)=noff+nr
1342 
1343  readbufferpointer(ioffp+nbuf)=noff ! pointer to start of buffer
1344  readbufferdatai(noff )=noff+nr ! pointer to end of buffer
1345  readbufferdatai(noff-1)=ifd(kfile)+jrec ! global record number (available with LOOP2)
1346  readbufferdataf(noff )=real(kfile) ! file number
1347  readbufferdataf(noff-1)=wfd(kfile) ! weight
1348 
1349  IF ((noff+nr+2+ndimbuf >= ndata*(iact+1)).OR.(nbuf >= npointer)) exit files ! buffer full
1350  ELSE
1351  !$OMP ATOMIC
1352  skippedrecords=skippedrecords+1
1353  cycle records
1354  END IF
1355 
1356  END DO records
1357 
1358  readbufferinfo(1,ithr)=-jfile ! flag eof
1359  IF (kfd(1,jfile) == 0) THEN
1360  print *, 'PEREAD: file ', kfile, 'read the first time, found',jrec,' records'
1361  kfd(1,jfile)=-jrec
1362  END IF
1363  ! take next file
1364  !$OMP CRITICAL
1365  IF (ifile < nfilb) THEN
1366  ifile=ifile+1
1367  jrec=0
1368  readbufferinfo(1,ithr)=ifile
1369  readbufferinfo(3,ithr)=jrec
1370  END IF
1371  !$OMP END CRITICAL
1372  jfile=readbufferinfo(1,ithr)
1373 
1374  END DO files
1375  !$OMP END PARALLEL
1376  ! compress pointers
1377  nrd=readbufferinfo(4,1) ! buffers from 1 .thread
1378  DO k=2,nthr
1379  iact =readbufferinfo(2,k)
1380  ioffp=iact*npointer
1381  nbuf=readbufferinfo(4,k)
1382  DO l=1,nbuf
1383  readbufferpointer(nrd+l)=readbufferpointer(ioffp+l)
1384  END DO
1385  nrd=nrd+nbuf
1386  END DO
1387 
1388  more=0
1389  DO k=1,nthr
1390  jfile=readbufferinfo(1,k)
1391  IF (jfile > 0) THEN ! no eof yet
1392  readbufferinfo(2,k)=more
1393  more=more+1
1394  ELSE
1395  ! no more files, thread retires
1396  readbufferinfo(1,k)=0
1397  readbufferinfo(2,k)=-1
1398  readbufferinfo(3,k)=0
1399  END IF
1400  END DO
1401  ! record limit ?
1402  IF (mxrec > 0.AND.(ntot+nrd) >= mxrec) THEN
1403  nrd=mxrec-ntot
1404  more=-1
1405  DO k=1,nthr
1406  jfile=readbufferinfo(1,k)
1407  IF (jfile > 0) THEN ! rewind files
1408  nrc=readbufferinfo(3,k)
1409  IF (kfd(1,jfile) == 0) kfd(1,jfile)=-nrc
1410  IF (kfile <= nfilf) THEN
1411  lun=kfile+10
1412  rewind lun
1413  ELSE
1414  lun=kfile-nfilf
1415 #ifdef READ_C_FILES
1416  CALL resetc(lun)
1417 #endif
1418  END IF
1419  END IF
1420  END DO
1421  END IF
1422 
1423  ntot=ntot+nrd
1424  nrec=ntot
1425  numreadbuffer=nrd
1426 
1427  sumrecords=sumrecords+nrd
1428  minrecordsinblock=min(minrecordsinblock,nrd)
1429  maxrecordsinblock=max(maxrecordsinblock,nrd)
1430 
1431  DO WHILE (nloopn == 0.AND.ntot >= nrpr)
1432  WRITE(*,*) ' Record ',nrpr
1433  IF (nrpr < 100000) THEN
1434  nrpr=nrpr*10
1435  ELSE
1436  nrpr=nrpr+100000
1437  END IF
1438  END DO
1439 
1440  IF (ncache > 0.AND.nloopn <= 1.AND. npri < mpri.AND.mprint > 1) THEN
1441  npri=npri+1
1442  IF (npri == 1) WRITE(*,100)
1443  WRITE(*,101) nrec, nrd, more ,ifile
1444 100 format(/' PeRead records active file' &
1445  /' total block threads number')
1446 101 format(' PeRead',4i10)
1447  END IF
1448 
1449  IF (more <= 0) THEN
1450  ifile=0
1451  IF (floop) THEN
1452  ! check for file weights
1453  ds0=0.0d0
1454  ds1=0.0d0
1455  ds2=0.0d0
1456  maxrecordsize=0
1457  maxrecordfile=0
1458  DO k=1,nfilb
1459  IF (xfd(k) > maxrecordsize) THEN
1460  maxrecordsize=xfd(k)
1461  maxrecordfile=k
1462  END IF
1463  dw=dfloat(-kfd(1,k))
1464  IF (wfd(k) /= 1.0) nfilw=nfilw+1
1465  ds0=ds0+dw
1466  ds1=ds1+dw*dble(wfd(k))
1467  ds2=ds2+dw*dble(wfd(k)**2)
1468  END DO
1469  print *, 'PEREAD: file ', maxrecordfile, 'with max record size ', maxrecordsize
1470  IF (nfilw > 0.AND.ds0 > 0.0d0) THEN
1471  ds1=ds1/ds0
1472  ds2=ds2/ds0-ds1*ds1
1473  DO lun=6,lunlog,2
1474  WRITE(lun,177) nfilw,sngl(ds1),sngl(ds2)
1475 177 format(/' !!!!!',i4,' weighted binary files', &
1476  /' !!!!! mean, variance of weights =',2g12.4)
1477  END DO
1478  END IF
1479  ! integrate record numbers
1480  DO k=2,nfilb
1481  ifd(k)=ifd(k-1)-kfd(1,k-1)
1482  END DO
1483  ! sort
1484  IF (nthr > 1) CALL sort2k(kfd,nfilb)
1485  IF (skippedrecords > 0) THEN
1486  print *, 'PEREAD skipped records: ', skippedrecords
1487  ndimbuf=maxrecordsize/2 ! adjust buffer size
1488  END IF
1489  END IF
1490  lprint=.false.
1491  floop=.false.
1492  IF (ncache > 0.AND.nloopn <= 1.AND.mprint > 0) &
1493  WRITE(*,179) numblocks, sumrecords, minrecordsinblock, maxrecordsinblock
1494 179 format(/' Read cache usage (#blocks, #records, ', &
1495  'min,max records/block'/17x,4i10)
1496  END IF
1497  return
1498 
1499 END SUBROUTINE peread
1500 
1508 SUBROUTINE peprep(mode)
1509  USE mpmod
1510 
1511  IMPLICIT NONE
1512 
1513  INTEGER, INTENT(IN) :: mode
1514 
1515  INTEGER :: ibuf
1516  INTEGER :: ichunk
1517  INTEGER :: iproc
1518  INTEGER :: isfrst
1519  INTEGER :: islast
1520  INTEGER :: ist
1521  INTEGER :: j
1522  INTEGER :: ja
1523  INTEGER :: jb
1524  INTEGER :: jsp
1525  INTEGER :: nst
1526  INTEGER :: inone
1527  !$ INTEGER :: OMP_GET_THREAD_NUM
1528 
1529 
1530  isfrst(ibuf)=readbufferpointer(ibuf)+1
1531  islast(ibuf)=readbufferdatai(readbufferpointer(ibuf))
1532 
1533  IF (mode > 0) THEN
1534  ichunk=min((numreadbuffer+mthrd-1)/mthrd/32+1,256)
1535  ! parallelize record loop
1536  !$OMP PARALLEL DO &
1537  !$OMP DEFAULT(PRIVATE) &
1538  !$OMP SHARED(numReadBuffer,readBufferPointer,readBufferDataI,readBufferDataF,ICHUNK) &
1539  !$OMP SCHEDULE(DYNAMIC,ICHUNK)
1540  DO ibuf=1,numreadbuffer ! buffer for current record
1541  iproc=0
1542  !$ IPROC=OMP_GET_THREAD_NUM() ! thread number
1543  ist=isfrst(ibuf)
1544  nst=islast(ibuf)
1545  DO ! loop over measurements
1546  CALL isjajb(nst,ist,ja,jb,jsp)
1547  IF(jb == 0) exit
1548  DO j=1,ist-jb
1549  readbufferdatai(jb+j)=inone( readbufferdatai(jb+j) ) ! translate to index
1550  END DO
1551  END DO
1552  END DO
1553  !$OMP END PARALLEL DO
1554  END IF
1555 
1556  !$POMP INST BEGIN(peprep)
1557  IF (mode <= 0) THEN
1558  DO ibuf=1,numreadbuffer ! buffer for current record
1559  ist=isfrst(ibuf)
1560  nst=islast(ibuf)
1561  DO ! loop over measurements
1562  CALL isjajb(nst,ist,ja,jb,jsp)
1563  IF(jb == 0) exit
1564  DO j=1,ist-jb
1565  readbufferdatai(jb+j)=inone( readbufferdatai(jb+j) ) ! translate to index
1566  END DO
1567  END DO
1568  END DO
1569  END IF
1570  !$POMP INST END(peprep)
1571 
1572 END SUBROUTINE peprep
1573 
1596 SUBROUTINE isjajb(nst,is,ja,jb,jsp)
1597  USE mpmod
1598 
1599  IMPLICIT NONE
1600  REAL :: glder
1601  INTEGER :: i
1602  INTEGER :: inder
1603 
1604  INTEGER, INTENT(IN) :: nst
1605  INTEGER, INTENT(IN OUT) :: is
1606  INTEGER, INTENT(OUT) :: ja
1607  INTEGER, INTENT(OUT) :: jb
1608  INTEGER, INTENT(OUT) :: jsp
1609  SAVE
1610  ! ...
1611  inder(i)=readbufferdatai(i)
1612  glder(i)=readbufferdataf(i)
1613 
1614  jsp=0
1615  DO
1616  ja=0
1617  jb=0
1618  IF(is >= nst) return
1619  DO
1620  is=is+1
1621  IF(inder(is) == 0) exit
1622  END DO
1623  ja=is
1624  DO
1625  is=is+1
1626  IF(inder(is) == 0) exit
1627  END DO
1628  jb=is
1629  DO WHILE(inder(is+1) /= 0.AND.is < nst)
1630  is=is+1
1631  END DO
1632  IF(ja+1 /= jb.OR.glder(jb) >= 0.0) exit
1633  ! special data
1634  jsp=jb ! pointer to special data
1635  is=is+ifix(-glder(jb)+0.5) ! skip NSP words
1636  END DO
1637 
1638 END SUBROUTINE isjajb
1639 
1640 
1641 !***********************************************************************
1642 ! LOOPN ...
1648 
1649 SUBROUTINE loopn
1650  USE mpmod
1651 
1652  IMPLICIT NONE
1653  REAL :: dsum
1654  REAL :: elmt
1655  REAL :: factrj
1656  REAL :: factrk
1657  REAL :: glder
1658  REAL :: peakd
1659  REAL :: peaki
1660  REAL :: ratae
1661  REAL :: rhs
1662  REAL :: rloop
1663  REAL :: sgm
1664  REAL :: used
1665  REAL :: usei
1666  REAL :: weight
1667  INTEGER :: i
1668  INTEGER :: ia
1669  INTEGER :: ib
1670  INTEGER :: ibuf
1671  INTEGER :: inder
1672  INTEGER :: ioffb
1673  INTEGER :: ipr
1674  INTEGER :: isfrst
1675  INTEGER :: islast
1676  INTEGER :: itgbi
1677  INTEGER :: itgbij
1678  INTEGER :: itgbik
1679  INTEGER :: ivgb
1680  INTEGER :: ivgbij
1681  INTEGER :: ivgbik
1682  INTEGER :: j
1683  INTEGER :: k
1684  INTEGER :: lastit
1685  INTEGER :: lun
1686  INTEGER :: ncrit
1687  INTEGER :: ndfs
1688  INTEGER :: ngras
1689  INTEGER :: nparl
1690  INTEGER :: nr
1691  INTEGER :: nrej
1692  INTEGER:: inone
1693 
1694  DOUBLE PRECISION:: adder
1695  DOUBLE PRECISION::funref
1696  DOUBLE PRECISION::dchi2s
1697  DOUBLE PRECISION::sndf
1698  INTEGER(KIND=large):: ii
1699  SAVE
1700  ! ...
1701  isfrst(ibuf)=readbufferpointer(ibuf)+1
1702  islast(ibuf)=readbufferdatai(readbufferpointer(ibuf))
1703  inder(i)=readbufferdatai(i)
1704  glder(i)=readbufferdataf(i)
1705  ! ----- book and reset ---------------------------------------------
1706  IF(nloopn == 0) THEN ! first call
1707  lastit=-1
1708  iitera=0
1709  END IF
1710 
1711  nloopn=nloopn+1 ! increase loop counter
1712  ndfsum=0
1713  sumndf=0.0d0
1714  funref=0.0d0
1715 
1716  IF(nloopn == 1) THEN ! book histograms for 1. iteration
1717  CALL gmpdef(1,4,'Function value in iterations')
1718  IF (metsol == 3) THEN ! extend to GMRES, i.e. 4?
1719  CALL gmpdef(2,3,'Number of MINRES steps vs iteration nr')
1720  END IF
1721  CALL hmpdef( 5,0.0,0.0,'Number of degrees of freedom')
1722  CALL hmpdef(11,0.0,0.0,'Number of local parameters')
1723  CALL hmpdef(23,0.0,0.0, 'SQRT of diagonal elements without presigma')
1724  CALL hmpdef(24,0.0,0.0, 'Log10 of off-diagonal elements')
1725  CALL hmpdef(25,0.0,0.0, 'Relative individual pre-sigma')
1726  CALL hmpdef(26,0.0,0.0, 'Relative global pre-sigma')
1727  END IF
1728 
1729 
1730  CALL hmpdef(3,-prange,prange, & ! book
1731  'Normalized residuals of single (global) measurement')
1732  CALL hmpdef(12,-prange,prange, & ! book
1733  'Normalized residuals of single (local) measurement')
1734  CALL hmpdef(13,-prange,prange, & ! book
1735  'Pulls of single (global) measurement')
1736  CALL hmpdef(14,-prange,prange, & ! book
1737  'Pulls of single (local) measurement')
1738  CALL hmpdef(4,0.0,0.0,'Chi^2/Ndf after local fit')
1739  CALL gmpdef(4,5,'location, dispersion (res.) vs record nr')
1740  CALL gmpdef(5,5,'location, dispersion (pull) vs record nr')
1741 
1742  ! WRITE(*,*) 'LOOPN ', NLOOPN, ' executing ICALCM=', ICALCM
1743 
1744  ! reset
1745 
1746  globalvector=0.0d0 ! reset rhs vector IGVEC
1747  IF(icalcm == 1) THEN
1748  globalmatd=0.0d0
1749  globalmatf=0.
1750  IF (metsol >= 3) matprecond=0.0d0
1751  END IF
1752 
1753  IF(nloopn == 2) CALL hmpdef(6,0.0,0.0,'Down-weight fraction')
1754 
1755  newite=.false.
1756  IF(iterat /= lastit) THEN ! new iteration
1757  newite=.true.
1758  funref=fvalue
1759  IF(nloopn > 1) THEN
1760  nrej=nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3)
1761  ! CALL MEND
1762  IF(iterat == 1) THEN
1763  chicut=chirem
1764  ELSE IF(iterat >= 1) THEN
1765  chicut=sqrt(chicut)
1766  IF(chicut /= 0.0.AND.chicut < 1.5) chicut=1.0
1767  IF(chicut /= 0.0.AND.nrej == 0) chicut=1.0
1768  END IF
1769  END IF
1770  ! WRITE(*,111) ! header line
1771  END IF
1772 
1773  DO i=0,3
1774  nrejec(i)=0 ! reset reject counter
1775  END DO
1776  DO k=3,6
1777  writebufferheader(k)=0 ! cache usage
1778  writebufferheader(-k)=0
1779  END DO
1780  ! statistics per binary file
1781  DO i=1,nfilb
1782  jfd(i)=0
1783  cfd(i)=0.0
1784  dfd(i)=0
1785  END DO
1786  ! ----- read next data ----------------------------------------------
1787  DO
1788  CALL peread(nr) ! read records
1789  CALL peprep(1) ! prepare records
1790  ndfs =0
1791  sndf =0.0d0
1792  dchi2s=0.0d0
1793  CALL loopbf(nrejec,ndfs,sndf,dchi2s,nfiles,jfd,cfd,dfd)
1794  ndfsum=ndfsum+ndfs
1795  sumndf=sumndf+sndf
1796  CALL addsum(dchi2s)
1797  IF (nr <= 0) exit ! next block of events ?
1798  END DO
1799  ! sum up RHS (over threads) once (reduction in LOOPBF: summation for each block)
1800  ioffb=0
1801  DO ipr=2,mthrd
1802  ioffb=ioffb+lenglobalvec
1803  DO k=1,lenglobalvec
1804  globalvector(k)=globalvector(k)+globalvector(ioffb+k)
1805  END DO
1806  END DO
1807 
1808  IF (icalcm == 1) THEN
1809  ! PRINT *, ' cache/w ',(writeBufferHeader(-K),K=3,6),(writeBufferHeader(K),K=3,6)
1810  nparl=writebufferheader(3)
1811  ncrit=writebufferheader(4)
1812  used=float(writebufferheader(-5))/float(writebufferheader(-3))*0.1
1813  usei=float(writebufferheader(5))/float(writebufferheader(3))*0.1
1814  peakd=float(writebufferheader(-6))*0.1
1815  peaki=float(writebufferheader(6))*0.1
1816  WRITE(*,111) nparl,ncrit,usei,used,peaki,peakd
1817 111 format(' Write cache usage (#flush,#overrun,<levels>,', &
1818  'peak(levels))'/2i7,',',4(f6.1,'%'))
1819  END IF
1820 
1821  ! ----- after end-of-data add contributions from pre-sigma ---------
1822 
1823  IF(nloopn == 1) THEN
1824  ! plot diagonal elements
1825  elmt=0.0
1826  DO i=1,nvgb ! diagonal elements
1827  ii=0
1828  IF(matsto == 1) THEN
1829  ii=i
1830  ii=(ii*ii+ii)/2
1831  END IF
1832  IF(matsto == 2) ii=i
1833  IF(matsto == 3) ii=i
1834  IF(ii /= 0) THEN
1835  elmt=sngl(globalmatd(ii))
1836  IF(elmt > 0.0) CALL hmpent(23,1.0/sqrt(elmt))
1837  END IF
1838  END DO
1839  END IF
1840 
1841 
1842 
1843  ! add pre-sigma contributions to matrix diagonal
1844 
1845  ! WRITE(*,*) 'Adding to diagonal ICALCM IND6',ICALCM,IND6
1846 
1847  IF(icalcm == 1) THEN
1848  DO ivgb=1,nvgb ! add evtl. pre-sigma
1849  ! WRITE(*,*) 'Index ',IVGB,IVGB,QM(IND6+IVGB)
1850  IF(globalparpreweight(ivgb) /= 0.0) THEN
1851  IF(ivgb > 0) CALL mupdat(ivgb,ivgb,dble(globalparpreweight(ivgb)))
1852  END IF
1853  END DO
1854  END IF
1855 
1856  CALL hmpwrt(23)
1857  CALL hmpwrt(24)
1858  CALL hmpwrt(25)
1859  CALL hmpwrt(26)
1860 
1861 
1862  ! add regularization term to F and to rhs --------------------------
1863 
1864  ! WRITE(*,*) 'NREGUL ',NREGUL,NLOOPN
1865 
1866  IF(nregul /= 0) THEN ! add regularization term to F and to rhs
1867  DO ivgb=1,nvgb
1868  itgbi=globalparvartototal(ivgb) ! global parameter index
1869  globalvector(ivgb)=globalvector(ivgb) -globalparameter(itgbi)*globalparpreweight(ivgb)
1870  adder=globalparpreweight(ivgb)*globalparameter(itgbi)**2
1871  CALL addsum(adder)
1872  END DO
1873  END IF
1874 
1875 
1876  ! ----- add contributions from "measurement" -----------------------
1877 
1878 
1879  i=1
1880  DO WHILE (i <= lenmeasurements)
1881  rhs=listmeasurements(i )%value ! right hand side
1882  sgm=listmeasurements(i+1)%value ! sigma parameter
1883  i=i+2
1884  weight=0.0
1885  IF(sgm > 0.0) weight=1.0/sgm**2
1886 
1887  dsum=-rhs
1888 
1889  ! loop over label/factor pairs
1890  ia=i
1891  DO
1892  i=i+1
1893  IF(i > lenmeasurements) exit
1894  IF(listmeasurements(i)%label == 0) exit
1895  END DO
1896  ib=i-1
1897 
1898  DO j=ia,ib
1899  factrj=listmeasurements(j)%value
1900  itgbij=inone(listmeasurements(j)%label) ! total parameter index
1901  IF(itgbij /= 0) THEN
1902  dsum=dsum+factrj*sngl(globalparameter(itgbij)) ! residuum
1903  END IF
1904  ! add to vector
1905  ivgbij=0
1906  IF(itgbij /= 0) ivgbij=globalparlabelindex(2,itgbij) ! variable-parameter index
1907  IF(ivgbij > 0) THEN
1908  globalvector(ivgbij)=globalvector(ivgbij) -weight*dsum*factrj ! vector
1909  END IF
1910 
1911  IF(icalcm == 1.AND.ivgbij > 0) THEN
1912  DO k=ia,j
1913  factrk=listmeasurements(k)%value
1914  itgbik=inone(listmeasurements(k)%label) ! total parameter index
1915  ! add to matrix
1916  ivgbik=0
1917  IF(itgbik /= 0) ivgbik=globalparlabelindex(2,itgbik) ! variable-parameter index
1918  IF(ivgbij > 0.AND.ivgbik > 0) THEN !
1919  CALL mupdat(ivgbij,ivgbik,dble(weight*factrj*factrk))
1920  END IF
1921  END DO
1922  END IF
1923  END DO
1924 
1925  adder=weight*dsum**2
1926  CALL addsum(adder) ! accumulate chi-square
1927 
1928  END DO
1929 
1930  ! ----- printout ---------------------------------------------------
1931 
1932 
1933  CALL getsum(fvalue) ! get accurate sum (Chi^2)
1934  flines=0.5d0*fvalue ! Likelihood function value
1935  rloop=iterat+0.01*nloopn
1936  actfun=sngl(funref-fvalue)
1937  IF(nloopn == 1) actfun=0.0
1938  ngras=nint(angras)
1939  ratae=0.0 !!!
1940  IF(delfun /= 0.0) THEN
1941  ratae=min(99.9,actfun/delfun) !!!
1942  ratae=max(-99.9,ratae)
1943  END IF
1944 
1945  ! rejects ...
1946 
1947  nrej =nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3)
1948  IF(nloopn == 1) THEN
1949  IF(nrej /= 0) THEN
1950  WRITE(*,*) ' '
1951  WRITE(*,*) 'Data rejected in initial loop:'
1952  WRITE(*,*) ' ', &
1953  nrejec(0), ' (rank deficit/NaN) ',nrejec(1),' (Ndf=0) ', &
1954  nrejec(2), ' (huge) ',nrejec(3),' (large)'
1955  END IF
1956  END IF
1957  ! IF(NREJEC(1)+NREJEC(2)+NREJEC(3).NE.0) THEN
1958  ! WRITE(LUNLOG,*) 'Data rejected in initial loop:',NREJEC(1),
1959  ! + ' (Ndf=0) ',NREJEC(2),' (huge) ',NREJEC(3),' (large)'
1960  ! END IF
1961 
1962 
1963  IF(newite.AND.iterat == 2) THEN
1964  IF(nrecpr /= 0.OR.nrecp2 /= 0) nrecer=nrec3
1965  IF(nrecpr < 0) THEN
1966  nrecpr=nrec1
1967  END IF
1968  IF(nrecp2 < 0) THEN
1969  nrecp2=nrec2
1970  END IF
1971  END IF
1972 
1973  IF(nloopn <= 2) THEN
1974  IF(nhistp /= 0) THEN
1975  ! CALL HMPRNT(3) ! scaled residual of single measurement
1976  ! CALL HMPRNT(12) ! scaled residual of single measurement
1977  ! CALL HMPRNT(4) ! chi^2/Ndf
1978  END IF
1979  CALL hmpwrt(3)
1980  CALL hmpwrt(12)
1981  CALL hmpwrt(4)
1982  CALL gmpwrt(4) ! location, dispersion (res.) as a function of record nr
1983  IF (nloopn <= lfitnp) THEN
1984  CALL hmpwrt(13)
1985  CALL hmpwrt(14)
1986  CALL gmpwrt(5) ! location, dispersion (pull) as a function of record nr
1987  END IF
1988  END IF
1989  ! IF(NLOOPN.EQ.2.AND.NHISTP.NE.0) CALL HMPRNT(6)
1990  IF(nloopn == 2) CALL hmpwrt(6)
1991  IF(nloopn <= 1) THEN
1992  ! IF(NHISTP.NE.0) CALL HMPRNT(5) ! number of degrees of freedom
1993  ! IF(NHISTP.NE.0) CALL HMPRNT(11) ! Nlocal
1994  CALL hmpwrt(5)
1995  CALL hmpwrt(11)
1996  END IF
1997 
1998  ! local fit: band matrix structure !?
1999  IF (nloopn == 1.AND.nbndr > 0) THEN
2000  DO lun=6,8,2
2001  WRITE(lun,*) ' '
2002  WRITE(lun,*) ' === local fits have bordered band matrix structure ==='
2003  WRITE(lun,101) ' NBNDR',nbndr,'number of records'
2004  WRITE(lun,101) ' NBDRX',nbdrx,'max border size'
2005  WRITE(lun,101) ' NBNDX',nbndx,'max band width'
2006  END DO
2007  END IF
2008 
2009  lastit=iterat
2010 
2011 101 format(1x,a6,' =',i10,' = ',a)
2012 ! 101 FORMAT(' LOOPN',I6,' Function value',F22.8,10X,I6,' records')
2013 ! 102 FORMAT(' incl. constraint penalty',F22.8)
2014 ! 103 FORMAT(I13,3X,A,G12.4)
2015 END SUBROUTINE loopn ! loop with fits
2016 
2020 
2021 SUBROUTINE ploopa(lunp)
2022  USE mpmod
2023 
2024  IMPLICIT NONE
2025 
2026  INTEGER, INTENT(IN) :: lunp
2027  ! ..
2028  WRITE(lunp,*) ' '
2029  WRITE(lunp,101) ! header line
2030  WRITE(lunp,102) ! header line
2031 101 format(' it fc',' fcn_value dfcn_exp slpr costh iit st', &
2032  ' ls step cutf',1x,'rejects hmmsec FMS')
2033 102 format(' -- --',' ----------- -------- ---- ----- --- --', &
2034  ' -- ---- ----',1x,'------- ------ ---')
2035  return
2036 END SUBROUTINE ploopa ! title for iteration
2037 
2041 
2042 SUBROUTINE ploopb(lunp)
2043  USE mpmod
2044 
2045  IMPLICIT NONE
2046  INTEGER :: ma
2047  INTEGER :: minut
2048  INTEGER :: nfa
2049  INTEGER :: nhour
2050  INTEGER :: nrej
2051  INTEGER :: nsecnd
2052  REAL :: ratae
2053  REAL :: rstb
2054  REAL :: secnd
2055  REAL :: slopes(3)
2056  REAL :: steps(3)
2057  REAL, DIMENSION(2) :: ta
2058 
2059  INTEGER, INTENT(IN) :: lunp
2060 
2061  CHARACTER (LEN=4):: ccalcm(4)
2062  DATA ccalcm / ' end',' S', ' F ',' FMS' /
2063  SAVE
2064 
2065  nrej=nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) ! rejects
2066  IF(nrej > 9999999) nrej=9999999
2067  CALL etime(ta,rstb)
2068  deltim=rstb-rstart
2069  CALL sechms(deltim,nhour,minut,secnd) ! time
2070  nsecnd=nint(secnd)
2071  IF(iterat == 0) THEN
2072  WRITE(lunp,103) iterat,nloopn,fvalue, &
2073  chicut,nrej,nhour,minut,nsecnd,ccalcm(lcalcm)
2074  ELSE
2075  CALL ptlopt(nfa,ma,slopes,steps) ! slopes steps
2076  ratae=abs(slopes(2)/slopes(1))
2077  stepl=steps(2)
2078  WRITE(lunp,104) iterat,nloopn,fvalue,delfun,ratae,angras, &
2079  iitera,istopa,lsinfo,stepl, chicut,nrej,nhour,minut,nsecnd,ccalcm(lcalcm)
2080  END IF
2081 103 format(i3,i3,e12.5,38x,f5.1, 1x,i7, i2,i2,i3,a4)
2082 104 format(i3,i3,e12.5,1x,e8.2,f6.3,f6.3,i5,2i3,f6.2,f5.1, &
2083  1x,i7, i2,i2,i3,a4)
2084  return
2085 END SUBROUTINE ploopb ! iteration line
2086 
2090 
2091 SUBROUTINE ploopc(lunp)
2092  USE mpmod
2093 
2094  IMPLICIT NONE
2095  INTEGER :: ma
2096  INTEGER :: minut
2097  INTEGER :: nfa
2098  INTEGER :: nhour
2099  INTEGER :: nrej
2100  INTEGER :: nsecnd
2101  REAL :: ratae
2102  REAL :: rstb
2103  REAL :: secnd
2104  REAL :: slopes(3)
2105  REAL :: steps(3)
2106  REAL, DIMENSION(2) :: ta
2107 
2108  INTEGER, INTENT(IN) :: lunp
2109  CHARACTER (LEN=4):: ccalcm(4)
2110  DATA ccalcm / ' end',' S', ' F ',' FMS' /
2111  SAVE
2112 
2113  nrej=nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) ! rejects
2114  IF(nrej > 9999999) nrej=9999999
2115  CALL etime(ta,rstb)
2116  deltim=rstb-rstart
2117  CALL sechms(deltim,nhour,minut,secnd) ! time
2118  nsecnd=nint(secnd)
2119  CALL ptlopt(nfa,ma,slopes,steps) ! slopes steps
2120  ratae=abs(slopes(2)/slopes(1))
2121  stepl=steps(2)
2122  WRITE(lunp,105) nloopn,fvalue, ratae,lsinfo, &
2123  stepl,nrej,nhour,minut,nsecnd,ccalcm(lcalcm)
2124 105 format(3x,i3,e12.5,9x, f6.3,14x,i3,f6.2,6x, i7, i2,i2,i3,a4)
2125  return
2126 
2127 END SUBROUTINE ploopc ! sub-iteration line
2128 
2132 
2133 SUBROUTINE ploopd(lunp)
2134  USE mpmod
2135  IMPLICIT NONE
2136  INTEGER :: minut
2137  INTEGER :: nhour
2138  INTEGER :: nsecnd
2139  REAL :: rstb
2140  REAL :: secnd
2141  REAL, DIMENSION(2) :: ta
2142 
2143  INTEGER, INTENT(IN) :: lunp
2144  CHARACTER (LEN=4):: ccalcm(4)
2145  DATA ccalcm / ' end',' S', ' F ',' FMS' /
2146  SAVE
2147  CALL etime(ta,rstb)
2148  deltim=rstb-rstart
2149  CALL sechms(deltim,nhour,minut,secnd) ! time
2150  nsecnd=ifix(secnd+0.5)
2151 
2152  WRITE(lunp,106) nhour,minut,nsecnd,ccalcm(lcalcm)
2153 106 format(69x,i2,i2,i3,a4)
2154  return
2155 END SUBROUTINE ploopd
2156 
2158 SUBROUTINE explfc(lunit)
2159  IMPLICIT NONE
2160  INTEGER:: lunit
2161  WRITE(lunit,*) ' '
2162  WRITE(lunit,102) 'Explanation of iteration table'
2163  WRITE(lunit,102) '=============================='
2164  WRITE(lunit,101) 'it', &
2165  'iteration number. Global parameters are improved for it > 0.'
2166  WRITE(lunit,102) 'First function evaluation is called iteraton 0.'
2167  WRITE(lunit,101) 'fc', 'number of function evaluations.'
2168  WRITE(lunit,101) 'fcn_value', 'value of 2 x Likelihood function (LF).'
2169  WRITE(lunit,102) 'The final value is the chi^2 value of the fit and should'
2170  WRITE(lunit,102) 'be about equal to the NDF (see below).'
2171  WRITE(lunit,101) 'dfcn_exp', &
2172  'expected reduction of the value of the Likelihood function (LF)'
2173  WRITE(lunit,101) 'slpr', 'ratio of the actual slope to inital slope.'
2174  WRITE(lunit,101) 'costh', &
2175  'cosine of angle between search direction and -gradient'
2176  WRITE(lunit,101) 'iit', &
2177  'number of internal iterations in GMRES/MINRES algorithmus'
2178  WRITE(lunit,101) 'st', 'stop code of GMRES/MINRES algorithmus'
2179  WRITE(lunit,102) '< 0: rhs is very special, with beta2 = 0'
2180  WRITE(lunit,102) '= 0: rhs b = 0, i.e. the exact solution is x = 0'
2181  WRITE(lunit,102) '= 1 requested accuracy achieved, as determined by rtol'
2182  WRITE(lunit,102) '= 2 reasonable accuracy achieved, given eps'
2183  WRITE(lunit,102) '= 3 x has converged to an eigenvector'
2184  WRITE(lunit,102) '= 4 matrix ill-conditioned (Acond has exceeded 0.1/eps)'
2185  WRITE(lunit,102) '= 5 the iteration limit was reached'
2186  WRITE(lunit,102) '= 6 Matrix x vector does not define a symmetric matrix'
2187  WRITE(lunit,102) '= 7 Preconditioner does not define a symmetric matrix'
2188  WRITE(lunit,101) 'ls', 'line search info'
2189  WRITE(lunit,102) '< 0 recalculate function'
2190  WRITE(lunit,102) '= 0: N or STP lt 0 or step not descending'
2191  WRITE(lunit,102) '= 1: Linesearch convergence conditions reached'
2192  WRITE(lunit,102) '= 2: interval of uncertainty at lower limit'
2193  WRITE(lunit,102) '= 3: max nr of line search calls reached'
2194  WRITE(lunit,102) '= 4: step at the lower bound'
2195  WRITE(lunit,102) '= 5: step at the upper bound'
2196  WRITE(lunit,102) '= 6: rounding error limitation'
2197  WRITE(lunit,101) 'step', &
2198  'the factor for the Newton step during the line search. Usually'
2199  WRITE(lunit,102) &
2200  'a value of 1 gives a sufficient reduction of the LF. Oherwise'
2201  WRITE(lunit,102) 'other step values are tried.'
2202  WRITE(lunit,101) 'cutf', &
2203  'cut factor. Local fits are rejected, if their chi^2 value'
2204  WRITE(lunit,102) &
2205  'is larger than the 3-sigma chi^2 value times the cut factor.'
2206  WRITE(lunit,102) 'A cut factor of 1 is used finally, but initially a larger'
2207  WRITE(lunit,102) 'factor may be used. A value of 0.0 means no cut.'
2208  WRITE(lunit,101) 'rejects', 'total number of rejected local fits.'
2209  WRITE(lunit,101) 'hmmsec', 'the time in hours (h), minutes (mm) and seconds.'
2210  WRITE(lunit,101) 'FMS', 'calculation of Function value, Matrix, Solution.'
2211  WRITE(lunit,*) ' '
2212 
2213 101 format(a9,' = ',a)
2214 102 format(13x,a)
2215 END SUBROUTINE explfc
2216 
2224 
2225 SUBROUTINE mupdat(i,j,add) !
2226  USE mpmod
2227 
2228  IMPLICIT NONE
2229 
2230  INTEGER, INTENT(IN) :: i
2231  INTEGER, INTENT(IN) :: j
2232  DOUBLE PRECISION, INTENT(IN) :: add
2233 
2234  INTEGER(KIND=large):: ijadd
2235  INTEGER(KIND=large):: ia
2236  INTEGER(KIND=large):: ja
2237  INTEGER(KIND=large):: ij
2238  ! ...
2239  IF(i <= 0.OR.j <= 0) return
2240  ia=max(i,j) ! larger
2241  ja=min(i,j) ! smaller
2242  ij=0
2243  IF(matsto == 1) THEN ! full symmetric matrix
2244  ij=ja+(ia*ia-ia)/2 ! ISYM index
2245  globalmatd(ij)=globalmatd(ij)+add
2246  ELSE IF(matsto == 2) THEN ! sparse symmetric matrix
2247  ij=ijadd(i,j) ! inline code requires same time
2248  IF (ij == 0) return ! pair is suppressed
2249  IF (ij > 0) THEN
2250  globalmatd(ij)=globalmatd(ij)+add
2251  ELSE
2252  globalmatf(-ij)=globalmatf(-ij)+sngl(add)
2253  END IF
2254  END IF
2255  IF(metsol >= 3) THEN
2256  IF(mbandw > 0) THEN ! for Cholesky decomposition
2257  IF(ia <= nvgb) THEN ! variable global parameter
2258  ij=indprecond(ia)-ia+ja
2259  IF(ia > 1.AND.ij <= indprecond(ia-1)) ij=0
2260  IF(ij /= 0) matprecond(ij)=matprecond(ij)+add
2261  IF(ij < 0.OR.ij > size(matprecond)) stop 'mupdat: bad index'
2262  ELSE ! Lagrange multiplier
2263  ij=indprecond(nvgb)+(ia-nvgb-1)*nvgb+ja
2264  IF(ij /= 0) matprecond(ij)=matprecond(ij)+add
2265  END IF
2266  ELSE IF(mbandw == 0) THEN ! default preconditioner
2267  IF(ia <= nvgb) THEN ! variable global parameter
2268  IF(ja == ia) matprecond(ia)=matprecond(ia)+add ! diag
2269  ELSE ! Lagrange multiplier
2270  ij=nvgb+(ia-nvgb-1)*nvgb+ja
2271  IF(ij /= 0) matprecond(ij)=matprecond(ij)+add
2272  END IF
2273  END IF
2274  END IF
2275 END SUBROUTINE mupdat
2276 
2277 
2307 
2308 SUBROUTINE loopbf(nrej,ndfs,sndf,dchi2s, numfil,naccf,chi2f,ndff)
2309  USE mpmod
2310 
2311  IMPLICIT NONE
2312  REAL :: cauchy
2313  REAL :: chichi
2314  REAL :: chlimt
2315  REAL :: chndf
2316  REAL :: chuber
2317  REAL :: down
2318  REAL :: glder
2319  REAL :: pull
2320  REAL :: r1
2321  REAL :: r2
2322  REAL :: rec
2323  REAL :: rerr
2324  REAL :: resid
2325  REAL :: resing
2326  REAL :: resmax
2327  REAL :: rmeas
2328  REAL :: rmloc
2329  REAL :: suwt
2330  REAL :: used
2331  REAL :: wght
2332  REAL :: wrc
2333  REAL :: chindl
2334  INTEGER :: i
2335  INTEGER :: ia
2336  INTEGER :: ib
2337  INTEGER :: ibuf
2338  INTEGER :: ichunk
2339  INTEGER :: icmn
2340  INTEGER :: icost
2341  INTEGER :: id
2342  INTEGER :: idiag
2343  INTEGER :: iext
2344  INTEGER :: ij
2345  INTEGER :: ije
2346  INTEGER :: ijn
2347  INTEGER :: ijsym
2348  INTEGER :: ik
2349  INTEGER :: ike
2350  INTEGER :: im
2351  INTEGER :: imeas
2352  INTEGER :: in
2353  INTEGER :: inder
2354  INTEGER :: inv
2355  INTEGER :: ioffb
2356  INTEGER :: ioffc
2357  INTEGER :: ioffd
2358  INTEGER :: ioffe
2359  INTEGER :: ioffi
2360  INTEGER :: iprdbg
2361  INTEGER :: iproc
2362  INTEGER :: irow
2363  INTEGER :: isfrst
2364  INTEGER :: islast
2365  INTEGER :: ist
2366  INTEGER :: iter
2367  INTEGER :: itgbi
2368  INTEGER :: ivgbj
2369  INTEGER :: ivgbk
2370  INTEGER :: j
2371  INTEGER :: ja
2372  INTEGER :: jb
2373  INTEGER :: jk
2374  INTEGER :: jn
2375  INTEGER :: joffd
2376  INTEGER :: joffi
2377  INTEGER :: jproc
2378  INTEGER :: jsp
2379  INTEGER :: k
2380  INTEGER :: kbdr
2381  INTEGER :: kbdrx
2382  INTEGER :: kbnd
2383  INTEGER :: kfl
2384  INTEGER :: kx
2385  INTEGER :: mbdr
2386  INTEGER :: mbnd
2387  INTEGER :: nalc
2388  INTEGER :: nalg
2389  INTEGER :: nan
2390  INTEGER :: nb
2391  INTEGER :: ndf
2392  INTEGER :: ndfsd
2393  INTEGER :: ndown
2394  INTEGER :: neq
2395  INTEGER :: nfred
2396  INTEGER :: nfrei
2397  INTEGER :: ngg
2398  INTEGER :: nprdbg
2399  INTEGER :: nrank
2400  INTEGER :: nrc
2401  INTEGER :: nst
2402  INTEGER :: nter
2403  INTEGER :: nweig
2404 
2405  INTEGER, INTENT(IN OUT) :: nrej(0:3)
2406  INTEGER, INTENT(IN OUT) :: ndfs
2407  DOUBLE PRECISION, INTENT(IN OUT) :: sndf
2408  DOUBLE PRECISION, INTENT(IN OUT) :: dchi2s
2409  INTEGER, INTENT(IN) :: numfil
2410  INTEGER, INTENT(IN OUT) :: naccf(numfil)
2411  REAL, INTENT(IN OUT) :: chi2f(numfil)
2412  INTEGER, INTENT(IN OUT) :: ndff(numfil)
2413 
2414  DOUBLE PRECISION::dch2sd
2415  DOUBLE PRECISION:: dchi2
2416  DOUBLE PRECISION::dvar
2417  DOUBLE PRECISION:: dw1
2418  DOUBLE PRECISION::dw2
2419  DOUBLE PRECISION::summ
2420 
2421  !$ INTEGER OMP_GET_THREAD_NUM
2422 
2423  LOGICAL:: lprnt
2424  LOGICAL::lhist
2425  CHARACTER (LEN=3):: chast
2426  DATA chuber/1.345/ ! constant for Huber down-weighting
2427  DATA cauchy/2.3849/ ! constant for Cauchy down-weighting
2428  SAVE chuber,cauchy
2429  ! ...
2430  ijsym(i,j)=min(i,j)+(max(i,j)*max(i,j)-max(i,j))/2
2431  isfrst(ibuf)=readbufferpointer(ibuf)+1
2432  islast(ibuf)=readbufferdatai(readbufferpointer(ibuf))
2433  inder(i)=readbufferdatai(i)
2434  glder(i)=readbufferdataf(i)
2435 
2436  ichunk=min((numreadbuffer+mthrd-1)/mthrd/32+1,256)
2437  ! reset header, 3 words per thread:
2438  ! number of entries, offset to data, indices
2439  writebufferinfo=0
2440  writebufferdata=0.
2441  ! IF (ICALCM.EQ.1) print *, ' ICHUNK ', MQ(IDOT3),MTHRD,ICHUNK
2442  nprdbg=0
2443  iprdbg=-1
2444 
2445  ! parallelize record loop
2446  ! private copy of NDFS,.. for each thread, combined at end, init with 0.
2447  !$OMP PARALLEL DO &
2448  !$OMP DEFAULT(PRIVATE) &
2449  !$OMP SHARED(numReadBuffer,readBufferPointer,readBufferDataI, &
2450  !$OMP readBufferDataF,writeBufferHeader,writeBufferInfo, &
2451  !$OMP writeBufferData,writeBufferIndices,writeBufferUpdates,globalVector, &
2452  !$OMP globalParameter,globalParLabelIndex,globalIndexUsage,backIndexUsage, &
2453  !$OMP NAGB,NVGB,NAGBN,ICALCM,ICHUNK,NLOOPN,NRECER,NPRDBG,IPRDBG, &
2454  !$OMP NEWITE,CHICUT,LHUBER,CHUBER,ITERAT,NRECPR,MTHRD, &
2455  !$OMP DWCUT,CHHUGE,NRECP2,CAUCHY,LFITNP,LFITBB) &
2456  !$OMP REDUCTION(+:NDFS,SNDF,DCHI2S,NREJ,NBNDR,NACCF,CHI2F,NDFF) &
2457  !$OMP REDUCTION(MAX:NBNDX,NBDRX) &
2458  !$OMP REDUCTION(MIN:NREC3) &
2459  !$OMP SCHEDULE(DYNAMIC,ICHUNK)
2460  DO ibuf=1,numreadbuffer ! buffer for current record
2461  nrc=readbufferdatai(isfrst(ibuf)-2) ! record
2462  kfl=nint(readbufferdataf(isfrst(ibuf)-1)) ! file
2463  wrc=readbufferdataf(isfrst(ibuf)-2) ! weight
2464  dw1=dble(wrc)
2465  dw2=dsqrt(dw1)
2466 
2467  iproc=0
2468  !$ IPROC=OMP_GET_THREAD_NUM() ! thread number
2469  ioffb=nagb*iproc ! offset 'f'.
2470  ioffc=nagbn*iproc ! offset 'c'.
2471  ioffe=nvgb*iproc ! offset 'e'
2472  ioffd=writebufferheader(-1)*iproc+writebufferinfo(2,iproc+1) ! offset data
2473  ioffi=writebufferheader(1)*iproc+writebufferinfo(3,iproc+1)+2 ! offset indices
2474  ! ----- reset ------------------------------------------------------
2475  lprnt=.false.
2476  lhist=(iproc == 0)
2477  rec=nrc ! floating point value
2478  IF(nloopn == 1.AND.mod(nrc,100000) == 0) THEN
2479  WRITE(*,*) 'Record',nrc,' ... still reading'
2480  END IF
2481 
2482  ! printout/debug only for one thread at a time
2483 
2484 
2485  ! flag for record printout -----------------------------------------
2486 
2487  lprnt=.false.
2488  IF(newite.AND.(iterat == 1.OR.iterat == 3)) THEN
2489  IF(nrc == nrecpr) lprnt=.true.
2490  IF(nrc == nrecp2) lprnt=.true.
2491  IF(nrc == nrecer) lprnt=.true.
2492  END IF
2493  IF (lprnt)THEN
2494  !$OMP ATOMIC
2495  nprdbg=nprdbg+1 ! number of threads with debug
2496  IF (nprdbg == 1) iprdbg=iproc ! first thread with debug
2497  IF (iproc /= iprdbg) lprnt=.false.
2498  ! print *, ' LPRNT ', NRC, NPRDBG, IPRDBG, IPROC, LPRNT
2499  END IF
2500  IF(lprnt) THEN
2501  WRITE(1,*) ' '
2502  WRITE(1,*) '------------------ Loop',nloopn, &
2503  ': Printout for record',nrc,iproc
2504  WRITE(1,*) ' '
2505  END IF
2506 
2507  ! ----- print data -------------------------------------------------
2508 
2509  IF(lprnt) THEN
2510  imeas=0 ! local derivatives
2511  ist=isfrst(ibuf)
2512  nst=islast(ibuf)
2513  DO ! loop over measurements
2514  CALL isjajb(nst,ist,ja,jb,jsp)
2515  IF(ja == 0) exit
2516  IF(imeas == 0) WRITE(1,1121)
2517  imeas=imeas+1
2518  WRITE(1,1122) imeas,glder(ja),glder(jb), &
2519  (inder(ja+j),glder(ja+j),j=1,jb-ja-1)
2520  END DO
2521 1121 format(/'Measured value and local derivatives'/ &
2522  ' i measured std_dev index...derivative ...')
2523 1122 format(i3,2g12.4,3(i3,g12.4)/(27x,3(i3,g12.4)))
2524 
2525  imeas=0 ! global derivatives
2526  ist=isfrst(ibuf)
2527  nst=islast(ibuf)
2528  DO ! loop over measurements
2529  CALL isjajb(nst,ist,ja,jb,jsp)
2530  IF(ja == 0) exit
2531  IF(imeas == 0) WRITE(1,1123)
2532  imeas=imeas+1
2533  IF (jb < ist) THEN
2534  IF(ist-jb > 2) THEN
2535  WRITE(1,1124) imeas,(globalparlabelindex(1,inder(jb+j)),inder(jb+j), &
2536  globalparlabelindex(2,inder(jb+j)),glder(jb+j),j=1,ist-jb)
2537  ELSE
2538  WRITE(1,1125) imeas,(globalparlabelindex(1,inder(jb+j)),inder(jb+j), &
2539  globalparlabelindex(2,inder(jb+j)),glder(jb+j),j=1,ist-jb)
2540  END IF
2541  END IF
2542  END DO
2543 1123 format(/'Global derivatives'/ &
2544  ' i label gindex vindex derivative ...')
2545 1124 format(i3,2(i9,i7,i7,g12.4)/(3x,2(i9,i7,i7,g12.4)))
2546 1125 format(i3,2(i9,i7,i7,g12.4))
2547  END IF
2548 
2549  ! ----- first loop -------------------------------------------------
2550  ! ------ prepare local fit ------
2551  ! count local and global derivates
2552  ! subtract actual alignment parameters from the measured data
2553 
2554  IF(lprnt) THEN
2555  WRITE(1,*) ' '
2556  WRITE(1,*) 'Data corrections using values of global parameters'
2557  WRITE(1,*) '=================================================='
2558  WRITE(1,101)
2559  END IF
2560  nalg=0 ! count number of global derivatives
2561  nalc=0 ! count number of local derivatives
2562  imeas=0
2563  ist=isfrst(ibuf)
2564  nst=islast(ibuf)
2565  DO ! loop over measurements
2566  CALL isjajb(nst,ist,ja,jb,jsp)
2567  IF(ja == 0) exit
2568  rmeas=glder(ja) ! data
2569  ! subtract global ... from measured value
2570  DO j=1,ist-jb ! global parameter loop
2571  itgbi=inder(jb+j) ! global parameter label
2572  rmeas=rmeas-glder(jb+j)*sngl(globalparameter(itgbi)) ! subtract !!! reversed
2573  IF (icalcm == 1) THEN
2574  ij=globalparlabelindex(2,itgbi) ! index of variable global parameter
2575  IF(ij > 0) THEN
2576  ijn=backindexusage(ioffe+ij) ! get index of index
2577  IF(ijn == 0) THEN ! not yet included
2578  nalg=nalg+1 ! count
2579  globalindexusage(ioffc+nalg)=ij ! store global index
2580  backindexusage(ioffe+ij)=nalg ! store back index
2581  END IF
2582  END IF
2583  END IF
2584  END DO
2585  IF(lprnt) THEN
2586  imeas=imeas+1
2587  IF (jb < ist) WRITE(1,102) imeas,glder(ja),rmeas,glder(jb)
2588  END IF
2589  readbufferdataf(ja)=rmeas ! global contribution subtracted
2590  DO j=1,jb-ja-1 ! local parameter loop
2591  ij=inder(ja+j)
2592  nalc=max(nalc,ij) ! number of local parameters
2593  END DO
2594  END DO
2595 101 format(' index measvalue corrvalue sigma')
2596 102 format(i6,2x,2g12.4,' +-',g12.4)
2597 
2598  IF(nalc <= 0) go to 90
2599 
2600  ngg=(nalg*nalg+nalg)/2
2601  IF (icalcm == 1) THEN
2602  DO k=1,nalg*nalc
2603  localglobalmatrix(k)=0.0d0 ! reset global-local matrix
2604  END DO
2605  writebufferindices(ioffi-1)=nrc ! index header:
2606  writebufferindices(ioffi )=nalg ! event number, number of global par
2607  CALL sort1k(globalindexusage(ioffc+1),nalg) ! sort global par.
2608  DO k=1,nalg
2609  iext=globalindexusage(ioffc+k)
2610  writebufferindices(ioffi+k)=iext ! global par indices
2611  backindexusage(ioffe+iext)=k ! update back index
2612  END DO
2613  DO k=1,ngg
2614  writebufferupdates(ioffd+k)=0.0d0 ! reset global-global matrix
2615  END DO
2616  END IF
2617  ! ----- iteration start and check ---------------------------------
2618 
2619  nter=1 ! first loop without down-weighting
2620  IF(nloopn /= 1.AND.lhuber /= 0) nter=lhuber
2621 
2622  ! check matrix for bordered band structure (MBDR+MBND+1 <= NALC)
2623  mbnd=-1
2624  mbdr=nalc
2625  DO i=1, nalc
2626  ibandh(i)=0
2627  END DO
2628  irow=1
2629  idiag=1
2630  ndfsd=0
2631 
2632  iter=0
2633  resmax=0.0
2634  DO WHILE(iter < nter) ! outlier suppresssion iteration loop
2635  iter=iter+1
2636  resmax=0.0
2637  IF(lprnt) THEN
2638  WRITE(1,*) ' '
2639  WRITE(1,*) 'Outlier-suppression iteration',iter,' of',nter
2640  WRITE(1,*) '=========================================='
2641  WRITE(1,*) ' '
2642  imeas=0
2643  END IF
2644 
2645  ! ----- second loop ------------------------------------------------
2646  ! accumulate normal equations for local fit and determine solution
2647  DO i=1,nalc
2648  blvec(i)=0.0d0 ! reset vector
2649  END DO
2650  DO i=1,(nalc*nalc+nalc)/2 ! GF: FIXME - not really, local parameter number...
2651  clmat(i)=0.0d0 ! (p)reset matrix
2652  END DO
2653  neq=0
2654  ndown=0
2655  nweig=0
2656  ist=isfrst(ibuf)
2657  nst=islast(ibuf)
2658  DO ! loop over measurements
2659  CALL isjajb(nst,ist,ja,jb,jsp)
2660  IF(ja == 0) exit
2661  rmeas=glder(ja) ! data
2662  rerr =glder(jb) ! ... and the error
2663  wght =1.0/rerr**2 ! weight from error
2664  neq=neq+1 ! count equation
2665  nweig=nweig+1
2666  resid=rmeas-localcorrections(neq) ! subtract previous fit
2667  IF(nloopn /= 1.AND.iter /= 1.AND.lhuber /= 0) THEN
2668  IF(iter <= 3) THEN
2669  IF(abs(resid) > chuber*rerr) THEN ! down-weighting
2670  wght=wght*chuber*rerr/abs(resid)
2671  ndown=ndown+1
2672  END IF
2673  ELSE ! Cauchy
2674  wght=wght/(1.0+(resid/rerr/cauchy)**2)
2675  END IF
2676  END IF
2677 
2678  IF(lprnt.AND.iter /= 1.AND.nter /= 1) THEN
2679  chast=' '
2680  IF(abs(resid) > chuber*rerr) chast='* '
2681  IF(abs(resid) > 3.0*rerr) chast='** '
2682  IF(abs(resid) > 6.0*rerr) chast='***'
2683  IF(imeas == 0) WRITE(1,*) 'Second loop: accumulate'
2684  IF(imeas == 0) WRITE(1,103)
2685  imeas=imeas+1
2686  down=1.0/sqrt(wght)
2687  r1=resid/rerr
2688  r2=resid/down
2689  WRITE(1,104) imeas,rmeas,resid,rerr,r1,chast,r2
2690  END IF
2691 103 format(' index corrvalue residuum sigma', &
2692  ' nresid cnresid')
2693 104 format(i6,2x,2g12.4,' +-',g12.4,f7.2,1x,a3,f8.2)
2694 
2695  DO j=1,jb-ja-1 ! normal equations, local parameter loop
2696  ij=inder(ja+j) ! local parameter index J
2697  blvec(ij)=blvec(ij)+dble(wght)*dble(rmeas)*dble(glder(ja+j))
2698  DO k=1,j
2699  ik=inder(ja+k) ! local parameter index K
2700  jk=ijsym(ij,ik) ! index in symmetric matrix
2701  clmat(jk)=clmat(jk) & ! force double precision
2702  +dble(wght)*dble(glder(ja+j))*dble(glder(ja+k))
2703  ! check for band matrix substructure
2704  IF (iter == 1) THEN
2705  id=iabs(ij-ik)+1
2706  im=min(ij,ik)
2707  ibandh(id)=max(ibandh(id),im)
2708  END IF
2709  END DO
2710  END DO
2711  END DO
2712  ! for non trivial fits check for bordered band matrix structure
2713  IF (iter == 1.AND.nalc > 5.AND.lfitbb > 0) THEN
2714  kx=-1
2715  kbdr=0
2716  kbdrx=0
2717  icmn=nalc**3 ! cost (*6) should improve by at least factor 2
2718  DO k=nalc,2,-1
2719  kbnd=k-2
2720  kbdr=max(kbdr,ibandh(k))
2721  icost=6*(nalc-kbdr)*(kbnd+kbdr+1)**2+2*kbdr**3
2722  IF (icost < icmn) THEN
2723  icmn=icost
2724  kx=k
2725  kbdrx=kbdr
2726  END IF
2727  END DO
2728  IF (kx > 0) THEN
2729  mbnd=kx-2
2730  mbdr=kbdrx
2731  END IF
2732  END IF
2733 
2734  IF (mbnd >= 0) THEN
2735  ! fast solution for border banded matrix (inverse for ICALCM>0)
2736  IF (nloopn == 1) THEN
2737  nbndr=nbndr+1
2738  nbdrx=max(nbdrx,mbdr)
2739  nbndx=max(nbndx,mbnd)
2740  END IF
2741 
2742  inv=0
2743  IF (nloopn <= lfitnp.AND.iter == 1) inv=1 ! band part of inverse (for pulls)
2744  IF (icalcm == 1.OR.lprnt) inv=2 ! complete inverse
2745  CALL sqmibb(clmat,blvec,nalc,mbdr,mbnd,inv,nrank, &
2746  vbnd,vbdr,aux,vbk,vzru,scdiag,scflag)
2747  ELSE
2748  ! full inversion and solution
2749  inv=2
2750  CALL sqminv(clmat,blvec,nalc,nrank,scdiag,scflag)
2751  END IF
2752  ! check for NaNs
2753  nan=0
2754  DO k=1, nalc
2755  IF ((.NOT.(blvec(k) <= 0.0d0)).AND. (.NOT.(blvec(k) > 0.0d0))) nan=nan+1
2756  END DO
2757 
2758  IF(lprnt) THEN
2759  WRITE(1,*) ' '
2760  WRITE(1,*) 'Parameter determination:',nalc,' parameters,', ' rank=',nrank
2761  WRITE(1,*) '-----------------------'
2762  IF(ndown /= 0) WRITE(1,*) ' ',ndown,' data down-weighted'
2763  WRITE(1,*) ' '
2764  END IF
2765 
2766  ! ----- third loop -------------------------------------------------
2767  ! calculate single residuals remaining after local fit and chi^2
2768 
2769  summ=0.0d0
2770  suwt=0.0
2771  neq=0
2772  imeas=0
2773  ist=isfrst(ibuf)
2774  nst=islast(ibuf)
2775  DO ! loop over measurements
2776  CALL isjajb(nst,ist,ja,jb,jsp)
2777  IF(ja == 0) exit
2778  rmeas=glder(ja) ! data (global contrib. subtracted)
2779  rerr =glder(jb) ! ... and the error
2780  wght =1.0/rerr**2 ! weight from error
2781  neq=neq+1 ! count equation
2782  rmloc=0.0 ! local fit result reset
2783  DO j=1,jb-ja-1 ! local parameter loop
2784  ij=inder(ja+j)
2785  rmloc=rmloc+glder(ja+j)*sngl(blvec(ij)) ! local fit result
2786  END DO
2787  localcorrections(neq)=rmloc ! save local fit result
2788  rmeas=rmeas-rmloc ! reduced to residual
2789 
2790  ! calculate pulls? (needs covariance matrix)
2791  IF(iter == 1.AND.inv > 0.AND.nloopn <= lfitnp) THEN
2792  dvar=0.0d0
2793  DO j=1,jb-ja-1
2794  ij=inder(ja+j)
2795  DO k=1,jb-ja-1
2796  ik=inder(ja+k)
2797  jk=ijsym(ij,ik)
2798  dvar=dvar+clmat(jk)*dble(glder(ja+j))*dble(glder(ja+k))
2799  END DO
2800  END DO
2801  ! some variance left to define a pull?
2802  IF (0.999999d0/dble(wght) > dvar) THEN
2803  pull=rmeas/sngl(dsqrt(1.0d0/dble(wght)-dvar))
2804  IF (lhist) THEN
2805  IF (jb < ist) THEN
2806  CALL hmpent(13,pull) ! histogram pull
2807  CALL gmpms(5,rec,pull)
2808  ELSE
2809  CALL hmpent(14,pull) ! histogram pull
2810  END IF
2811  END IF
2812  END IF
2813  END IF
2814 
2815  IF(iter == 1.AND.jb < ist.AND.lhist) &
2816  CALL gmpms(4,rec,rmeas/rerr) ! residual (with global deriv.)
2817 
2818  dchi2=wght*rmeas*rmeas
2819  ! DCHIT=DCHI2
2820  resid=rmeas
2821  IF(nloopn /= 1.AND.iter /= 1.AND.lhuber /= 0) THEN
2822  IF(iter <= 3) THEN
2823  IF(abs(resid) > chuber*rerr) THEN ! down-weighting
2824  wght=wght*chuber*rerr/abs(resid)
2825  dchi2=2.0*chuber*(abs(resid)/rerr-0.5*chuber)
2826  END IF
2827  ELSE
2828  wght=wght/(1.0+(resid/rerr/cauchy)**2)
2829  dchi2=log(1.0+(resid/rerr/cauchy)**2)*cauchy**2
2830  END IF
2831  END IF
2832 
2833  down=1.0/sqrt(wght)
2834 
2835  ! SUWT=SUWT+DCHI2/DCHIT
2836  suwt=suwt+rerr/down
2837  IF(lprnt) THEN
2838  chast=' '
2839  IF(abs(resid) > chuber*rerr) chast='* '
2840  IF(abs(resid) > 3.0*rerr) chast='** '
2841  IF(abs(resid) > 6.0*rerr) chast='***'
2842  IF(imeas == 0) WRITE(1,*) 'Third loop: single residuals'
2843  IF(imeas == 0) WRITE(1,105)
2844  imeas=imeas+1
2845  r1=resid/rerr
2846  r2=resid/down
2847  IF(resid < 0.0) r1=-r1
2848  IF(resid < 0.0) r2=-r2
2849  WRITE(1,106) imeas,glder(ja),rmeas,rerr,r1,chast,r2
2850  END IF
2851 105 format(' index corrvalue residuum sigma', &
2852  ' nresid cnresid')
2853 106 format(i6,2x,2g12.4,' +-',g12.4,f7.2,1x,a3,f8.2)
2854 
2855  IF(iter == nter) THEN
2856  readbufferdataf(ja)=rmeas ! store remaining residual
2857  resmax=max(resmax,abs(rmeas)/rerr)
2858  END IF
2859 
2860  IF(iter == 1.AND.lhist) THEN
2861  IF (jb < ist) THEN
2862  CALL hmpent( 3,rmeas/rerr) ! histogram norm residual
2863  ELSE
2864  CALL hmpent(12,rmeas/rerr) ! histogram norm residual
2865  END IF
2866  END IF
2867  summ=summ+dchi2 ! accumulate chi-square sum
2868  END DO
2869  ndf=neq-nrank
2870  ! external seed ?
2871  dch2sd=0.0d0
2872  resing=(float(nweig)-suwt)/float(nweig)
2873  IF (lhist) THEN
2874  IF(iter == 1) CALL hmpent( 5,float(ndf)) ! histogram Ndf
2875  IF(iter == 1) CALL hmpent(11,float(nalc)) ! histogram Nlocal
2876  IF(nloopn == 2.AND.iter == nter) CALL hmpent(6,resing)
2877  END IF
2878  IF(lprnt) THEN
2879  WRITE(1,*) ' '
2880  WRITE(1,*) 'Chi^2=',summ,' at',ndf,' degrees of freedom: ', &
2881  '3-sigma limit is',chindl(3,ndf)*float(ndf)
2882  WRITE(1,*) suwt,' is sum of factors, compared to',nweig, &
2883  ' Downweight fraction:',resing
2884  END IF
2885  IF(nrank /= nalc.OR.nan > 0) THEN
2886  nrej(0)=nrej(0)+1 ! count cases
2887  IF (nrec3 == maxi4) nrec3=nrc
2888  IF(lprnt) THEN
2889  WRITE(1,*) ' rank deficit/NaN ', nalc, nrank, nan
2890  WRITE(1,*) ' ---> rejected!'
2891  END IF
2892  go to 90
2893  END IF
2894  IF(ndf <= 0) THEN
2895  nrej(1)=nrej(1)+1 ! count cases
2896  IF(lprnt) THEN
2897  WRITE(1,*) ' ---> rejected!'
2898  END IF
2899  go to 90
2900  END IF
2901 
2902  chndf=sngl(summ/dfloat(ndf))
2903 
2904  IF(iter == 1.AND.lhist) CALL hmpent(4,chndf) ! histogram chi^2/Ndf
2905  END DO ! outlier iteration loop
2906 
2907  ndfs=ndfs+ndf ! (local) sum of Ndf
2908  sndf=sndf+dfloat(ndf)*dw1 ! (local) weighted sum of Ndf
2909 
2910  ! ----- reject eventually ------------------------------------------
2911 
2912  IF(newite.AND.iterat == 2) THEN ! find record with largest Chi^2/Ndf
2913  IF(nrecp2 < 0.AND.chndf > writebufferdata(2,iproc+1)) THEN
2914  writebufferdata(2,iproc+1)=chndf
2915  writebufferinfo(7,iproc+1)=nrc
2916  END IF
2917  END IF
2918 
2919  chichi=chindl(3,ndf)*float(ndf)
2920  ! GF IF(SUMM.GT.50.0*CHICHI) THEN ! huge
2921  ! CHK CHICUT<0: NO cut (1st iteration)
2922  IF(chicut >= 0.0) THEN
2923  IF(summ > chhuge*chichi) THEN ! huge
2924  nrej(2)=nrej(2)+1 ! count cases with huge chi^2
2925  IF(lprnt) THEN
2926  WRITE(1,*) ' ---> rejected!'
2927  END IF
2928  go to 90
2929  END IF
2930 
2931  IF(chicut > 0.0) THEN
2932  chlimt=chicut*chichi
2933  ! WRITE(*,*) 'chi^2 ',SUMM,CHLIMT,CHICUT,CHINDL(3,NDF),NDF
2934  IF(summ > chlimt) THEN
2935  IF(lprnt) THEN
2936  WRITE(1,*) ' ---> rejected!'
2937  END IF
2938  ! add to FVALUE
2939  dchi2=chlimt ! total contribution limit
2940  dchi2s=dchi2s+dchi2*dw1 ! add total contribution
2941  nrej(3)=nrej(3)+1 ! count cases with large chi^2
2942  go to 90
2943  END IF
2944  END IF
2945  END IF
2946 
2947  IF(lhuber > 1.AND.dwcut /= 0.0.AND.resing > dwcut) THEN
2948  ! add to FVALUE
2949  dchi2=summ ! total contribution
2950  dchi2s=dchi2s+dchi2*dw1 ! add total contribution
2951  nrej(3)=nrej(3)+1 ! count cases with large chi^2
2952  ! WRITE(*,*) 'Downweight fraction cut ',RESING,DWCUT,SUMM
2953  IF(lprnt) THEN
2954  WRITE(1,*) ' ---> rejected!'
2955  END IF
2956  go to 90
2957  END IF
2958 
2959  IF(newite.AND.iterat == 2) THEN ! find record with largest residual
2960  IF(nrecpr < 0.AND.resmax > writebufferdata(1,iproc+1)) THEN
2961  writebufferdata(1,iproc+1)=resmax
2962  writebufferinfo(6,iproc+1)=nrc
2963  END IF
2964  END IF
2965  ! 'track quality' per binary file: accepted records
2966  naccf(kfl)=naccf(kfl)+1
2967  ndff(kfl) =ndff(kfl) +ndf
2968  chi2f(kfl)=chi2f(kfl)+chndf
2969 
2970  ! ----- fourth loop ------------------------------------------------
2971  ! update of global matrix and vector according to the "Millepede"
2972  ! principle, from the global/local information
2973 
2974  dchi2s=dchi2s+dch2sd
2975 
2976  ist=isfrst(ibuf)
2977  nst=islast(ibuf)
2978  DO ! loop over measurements
2979  CALL isjajb(nst,ist,ja,jb,jsp)
2980  IF(ja <= 0) exit
2981 
2982  rmeas=glder(ja) ! data residual
2983  rerr =glder(jb) ! ... and the error
2984  wght =1.0/rerr**2 ! weight from measurement error
2985  dchi2=wght*rmeas*rmeas ! least-square contribution
2986 
2987  IF(nloopn /= 1.AND.lhuber /= 0) THEN ! check residual
2988  resid=abs(rmeas)
2989  IF(resid > chuber*rerr) THEN
2990  wght=wght*chuber*rerr/resid ! down-weighting
2991  dchi2=2.0*chuber*(resid/rerr-0.5*chuber) ! modified contribution
2992  END IF
2993  END IF
2994  dchi2s=dchi2s+dchi2*dw1 ! add to total objective function
2995 
2996  ! global-global matrix contribution: add directly to gg-matrix
2997 
2998  DO j=1,ist-jb
2999  ivgbj=globalparlabelindex(2,inder(jb+j)) ! variable-parameter index
3000  IF(ivgbj > 0) THEN
3001  globalvector(ioffb+ivgbj)=globalvector(ioffb+ivgbj) &
3002  +dw1*wght*rmeas*glder(jb+j) ! vector !!! reverse
3003  IF(icalcm == 1) THEN
3004  ije=backindexusage(ioffe+ivgbj) ! get index of index, non-zero
3005  DO k=1,j
3006  ivgbk=globalparlabelindex(2,inder(jb+k))
3007  IF(ivgbk > 0) THEN
3008  ike=backindexusage(ioffe+ivgbk) ! get index of index, non-zero
3009  ia=max(ije,ike) ! larger
3010  ib=min(ije,ike) ! smaller
3011  ij=ib+(ia*ia-ia)/2
3012  writebufferupdates(ioffd+ij)=writebufferupdates(ioffd+ij) &
3013  -dw1*wght*glder(jb+j)*glder(jb+k)
3014  END IF
3015  END DO
3016  END IF
3017  END IF
3018  END DO
3019 
3020  ! normal equations - rectangular matrix for global/local pars
3021  ! global-local matrix contribution: accumulate rectangular matrix
3022  IF (icalcm /= 1) cycle
3023  DO j=1,ist-jb
3024  ivgbj=globalparlabelindex(2,inder(jb+j)) ! variable-parameter index
3025  IF(ivgbj > 0) THEN
3026  ije=backindexusage(ioffe+ivgbj) ! get index of index, non-zero
3027  DO k=1,jb-ja-1
3028  ik=inder(ja+k) ! local index
3029  jk=ik+(ije-1)*nalc ! matrix index
3030  localglobalmatrix(jk)=localglobalmatrix(jk)+dw2*wght*glder(jb+j)*glder(ja+k)
3031  END DO
3032  END IF
3033  END DO
3034  END DO
3035 
3036 
3037  ! ----- final matrix update ----------------------------------------
3038  ! update global matrices and vectors
3039  IF(icalcm /= 1) go to 90 ! matrix update
3040  ! (inverse local matrix) * (rectang. matrix) -> CORM
3041  ! T
3042  ! resulting symmetrix matrix = G * Gamma^{-1} * G
3043 
3044  CALL dbavat(clmat,localglobalmatrix,writebufferupdates(ioffd+1),nalc,-nalg)
3045 
3046  ! (rectang. matrix) * (local param vector) -> CORV
3047  ! resulting vector = G * q (q = local parameter)
3048  ! CALL DBGAX(DQ(IGLMA/2+1),BLVEC,DQ(ICORV/2+1),NALG,NALC) ! not done
3049  ! the vector update is not done, because after local fit it is zero!
3050 
3051  ! update cache status
3052  writebufferinfo(1,iproc+1)=writebufferinfo(1,iproc+1)+1
3053  writebufferinfo(2,iproc+1)=writebufferinfo(2,iproc+1)+ngg
3054  writebufferinfo(3,iproc+1)=writebufferinfo(3,iproc+1)+nalg+2
3055  ! check free space
3056  nfred=writebufferheader(-1)-writebufferinfo(2,iproc+1)-writebufferheader(-2)
3057  nfrei=writebufferheader(1)-writebufferinfo(3,iproc+1)-writebufferheader(2)
3058  IF (nfred < 0.OR.nfrei < 0) THEN ! need to flush
3059  nb=writebufferinfo(1,iproc+1)
3060  joffd=writebufferheader(-1)*iproc ! offset data
3061  joffi=writebufferheader(1)*iproc+2 ! offset indices
3062  used=float(writebufferinfo(2,iproc+1))/float(writebufferheader(-1))
3063  writebufferinfo(4,iproc+1)=writebufferinfo(4,iproc+1) +ifix(1000.0*used+0.5)
3064  used=float(writebufferinfo(3,iproc+1))/float(writebufferheader(1))
3065  writebufferinfo(5,iproc+1)=writebufferinfo(5,iproc+1) +ifix(1000.0*used+0.5)
3066  !$OMP CRITICAL
3067  writebufferheader(-4)=writebufferheader(-4)+1
3068  writebufferheader(4)=writebufferheader(4)+1
3069 
3070  DO ib=1,nb
3071  ijn=0
3072  DO in=1,writebufferindices(joffi)
3073  i=writebufferindices(joffi+in)
3074  ! DQ(IGVEC/2+I)=DQ(IGVEC/2+I)+DQ(ICORV/2+IN) ! not done: = zero
3075  DO jn=1,in
3076  ijn=ijn+1
3077  j=writebufferindices(joffi+jn)
3078  CALL mupdat(i,j,-writebufferupdates(joffd+ijn)) ! matrix update
3079  END DO
3080  END DO
3081  joffd=joffd+ijn
3082  joffi=joffi+writebufferindices(joffi)+2
3083  END DO
3084  !$OMP END CRITICAL
3085  ! reset counter, pointers
3086  DO k=1,3
3087  writebufferinfo(k,iproc+1)=0
3088  END DO
3089  END IF
3090 
3091 90 IF(lprnt) THEN
3092  WRITE(1,*) ' '
3093  WRITE(1,*) '------------------ End of printout for record',nrc
3094  WRITE(1,*) ' '
3095  END IF
3096 
3097  DO i=1,nalg ! reset global index array
3098  iext=globalindexusage(ioffc+i)
3099  backindexusage(ioffe+iext)=0
3100  END DO
3101 
3102  END DO
3103  !$OMP END PARALLEL DO
3104 
3105  IF (icalcm == 1) THEN
3106  ! flush remaining matrices
3107  DO k=1,mthrd ! update statistics
3108  writebufferheader(-3)=writebufferheader(-3)+1
3109  used=float(writebufferinfo(2,k))/float(writebufferheader(-1))
3110  writebufferinfo(4,k)=writebufferinfo(4,k)+ifix(1000.0*used+0.5)
3111  writebufferheader(-5)=writebufferheader(-5)+writebufferinfo(4,k)
3112  writebufferheader(-6)=max(writebufferheader(-6),writebufferinfo(4,k))
3113  writebufferinfo(4,k)=0
3114  writebufferheader(3)=writebufferheader(3)+1
3115  used=float(writebufferinfo(3,k))/float(writebufferheader(1))
3116  writebufferinfo(5,k)=writebufferinfo(5,k)+ifix(1000.0*used+0.5)
3117  writebufferheader(5)=writebufferheader(5)+writebufferinfo(5,k)
3118  writebufferheader(6)=max(writebufferheader(6),writebufferinfo(5,k))
3119  writebufferinfo(5,k)=0
3120  END DO
3121 
3122  !$OMP PARALLEL &
3123  !$OMP DEFAULT(PRIVATE) &
3124  !$OMP SHARED(writeBufferHeader,writeBufferInfo,writeBufferIndices,writeBufferUpdates,MTHRD)
3125  iproc=0
3126  !$ IPROC=OMP_GET_THREAD_NUM() ! thread number
3127  DO jproc=0,mthrd-1
3128  nb=writebufferinfo(1,jproc+1)
3129  ! print *, ' flush end ', JPROC, NRC, NB
3130  joffd=writebufferheader(-1)*jproc ! offset data
3131  joffi=writebufferheader(1)*jproc+2 ! offset indices
3132  DO ib=1,nb
3133  ! print *, ' buf end ', JPROC,IB,writeBufferIndices(JOFFI-1),writeBufferIndices(JOFFI)
3134  ijn=0
3135  DO in=1,writebufferindices(joffi)
3136  i=writebufferindices(joffi+in)
3137  !$ IF (MOD(I,MTHRD).EQ.IPROC) THEN
3138  DO jn=1,in
3139  ijn=ijn+1
3140  j=writebufferindices(joffi+jn)
3141  CALL mupdat(i,j,-writebufferupdates(joffd+ijn)) ! matrix update
3142  END DO
3143  !$ ELSE
3144  !$ IJN=IJN+IN
3145  !$ ENDIF
3146  END DO
3147  joffd=joffd+ijn
3148  joffi=joffi+writebufferindices(joffi)+2
3149  END DO
3150  END DO
3151  !$OMP END PARALLEL
3152  END IF
3153 
3154  IF(newite.AND.iterat == 2) THEN ! get worst records (for printrecord -1 -1)
3155  IF (nrecpr < 0) THEN
3156  DO k=1,mthrd
3157  IF (writebufferdata(1,k) > value1) THEN
3158  value1=writebufferdata(1,k)
3159  nrec1 =writebufferinfo(6,k)
3160  END IF
3161  END DO
3162  END IF
3163  IF (nrecp2 < 0) THEN
3164  DO k=1,mthrd
3165  IF (writebufferdata(2,k) > value2) THEN
3166  value2=writebufferdata(2,k)
3167  nrec2 =writebufferinfo(7,k)
3168  END IF
3169  END DO
3170  END IF
3171  END IF
3172 
3173 END SUBROUTINE loopbf
3174 
3175 
3176 
3177 
3178 !***********************************************************************
3179 
3191 SUBROUTINE prtglo
3192  USE mpmod
3193 
3194  IMPLICIT NONE
3195  REAL:: dpa
3196  REAL:: err
3197  REAL:: gcor
3198  INTEGER:: i
3199  INTEGER:: ie
3200  INTEGER:: iev
3201  INTEGER:: ij
3202  INTEGER:: imin
3203  INTEGER:: iprlim
3204  INTEGER:: isub
3205  INTEGER:: itgbi
3206  INTEGER:: itgbl
3207  INTEGER:: ivgbi
3208  INTEGER:: j
3209  INTEGER:: label
3210  INTEGER:: lup
3211  REAL:: par
3212 
3213  DOUBLE PRECISION:: diag
3214  DOUBLE PRECISION::gmati
3215  DOUBLE PRECISION::gcor2
3216  INTEGER:: labele(3)
3217  INTEGER(KIND=large):: ii
3218  REAL:: compnt(3)
3219  SAVE
3220  ! ...
3221 
3222  lup=09
3223  CALL mvopen(lup,'millepede.res')
3224 
3225  WRITE(*,*) ' '
3226  WRITE(*,*) ' Result of fit for global parameters'
3227  WRITE(*,*) ' ==================================='
3228  WRITE(*,*) ' '
3229 
3230  WRITE(*,101)
3231 
3232  WRITE(lup,*) 'Parameter ! first 3 elements per line are', &
3233  ' significant (if used as input)'
3234  iprlim=10
3235  DO itgbi=1,ntgb ! all parameter variables
3236  itgbl=globalparlabelindex(1,itgbi)
3237  ivgbi=globalparlabelindex(2,itgbi)
3238  par=sngl(globalparameter(itgbi)) ! initial value
3239  IF(ivgbi > 0) THEN
3240  dpa=par-globalparstart(itgbi) ! difference
3241  IF(metsol == 1.OR.metsol == 2) THEN
3242  ii=ivgbi
3243  ii=(ii*ii+ii)/2
3244  gmati=globalmatd(ii)
3245  err=sqrt(abs(sngl(gmati)))
3246  IF(gmati < 0.0d0) err=-err
3247  diag=workspaced(ivgbi)
3248  gcor=-1.0
3249  IF(gmati*diag > 0.0d0) THEN ! global correlation
3250  gcor2=1.0d0-1.0d0/(gmati*diag)
3251  IF(gcor2 >= 0.0d0.AND.gcor2 <= 1.0d0) gcor=sngl(dsqrt(gcor2))
3252  END IF
3253  END IF
3254  END IF
3255  IF(itgbi <= iprlim) THEN
3256  IF(ivgbi <= 0) THEN
3257  WRITE(* ,102) itgbl,par,globalparpresigma(itgbi)
3258  ELSE
3259  IF(metsol == 1.OR.metsol == 2) THEN
3260  IF (igcorr == 0) THEN
3261  WRITE(*,102) itgbl,par,globalparpresigma(itgbi),dpa,err
3262  ELSE
3263  WRITE(*,102) itgbl,par,globalparpresigma(itgbi),dpa,err,gcor
3264  END IF
3265  ELSE
3266  WRITE(*,102) itgbl,par,globalparpresigma(itgbi),dpa
3267  END IF
3268  END IF
3269  ELSE IF(itgbi == iprlim+1) THEN
3270  WRITE(* ,*) '... (further printout suppressed, but see log file)'
3271  END IF
3272 
3273  ! file output
3274  IF(ivgbi <= 0) THEN
3275  WRITE(lup,102) itgbl,par,globalparpresigma(itgbi)
3276  ELSE
3277  IF(metsol == 1.OR.metsol == 2) THEN
3278  IF (igcorr == 0) THEN
3279  WRITE(lup,102) itgbl,par,globalparpresigma(itgbi),dpa,err
3280  ELSE
3281  WRITE(lup,102) itgbl,par,globalparpresigma(itgbi),dpa,err,gcor
3282  END IF
3283  ELSE
3284  WRITE(lup,102) itgbl,par,globalparpresigma(itgbi),dpa
3285  END IF
3286  END IF
3287  END DO
3288  rewind lup
3289  CLOSE(unit=lup)
3290 
3291  IF(metsol == 2) THEN ! diagonalisation: write eigenvectors
3292  CALL mvopen(lup,'millepede.eve')
3293  imin=1
3294  DO i=nagb,1,-1
3295  IF(workspaceeigenvalues(i) > 0.0d0) THEN
3296  imin=i ! index of smallest pos. eigenvalue
3297  exit
3298  ENDIF
3299  END DO
3300  iev=0
3301 
3302  DO isub=0,min(15,imin-1)
3303  IF(isub < 10) THEN
3304  i=imin-isub
3305  ELSE
3306  i=isub-9
3307  END IF
3308 
3309  ! DO I=IMIN,MAX(1,IMIN-9),-1 ! backward loop, up to 10 vectors
3310  WRITE(*,*) 'Eigenvector ',i,' with eigenvalue',workspaceeigenvalues(i)
3311  WRITE(lup,*) 'Eigenvector ',i,' with eigenvalue',workspaceeigenvalues(i)
3312  DO j=1,nagb
3313  ij=j+(i-1)*nagb ! index with eigenvector array
3314  IF(j <= nvgb) THEN
3315  itgbi=globalparvartototal(j)
3316  label=globalparlabelindex(1,itgbi)
3317  ELSE
3318  label=nvgb-j ! label negative for constraints
3319  END IF
3320  iev=iev+1
3321  labele(iev)=label
3322  compnt(iev)=sngl(workspaceeigenvectors(ij)) ! component
3323  IF(iev == 3) THEN
3324  WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
3325  iev=0
3326  END IF
3327  END DO
3328  IF(iev /= 0) WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
3329  iev=0
3330  WRITE(lup,*) ' '
3331  END DO
3332 
3333  END IF
3334 
3335 101 format(1x,' label parameter presigma differ', &
3336  ' error'/ 1x,'-----------',4x,4('-------------'))
3337 102 format(i10,2x,4g14.5,f8.3)
3338 103 format(3(i11,f11.7,2x))
3339 END SUBROUTINE prtglo ! print final log file
3340 
3341 
3350 
3351 SUBROUTINE avprod(n,x,b)
3352  USE mpmod
3353 
3354  IMPLICIT NONE
3355  INTEGER :: i
3356  INTEGER :: iencdb
3357  INTEGER :: iencdm
3358  INTEGER :: iproc
3359  INTEGER :: ir
3360  INTEGER :: j
3361  INTEGER :: jc
3362  INTEGER :: jj
3363  INTEGER :: jn
3364 
3365  INTEGER, INTENT(IN) :: n
3366  DOUBLE PRECISION, INTENT(IN) :: x(n)
3367  DOUBLE PRECISION, INTENT(OUT) :: b(n)
3368  INTEGER(KIND=large) :: k
3369  INTEGER(KIND=large) :: kk
3370  INTEGER(KIND=large) :: kl
3371  INTEGER(KIND=large) :: ku
3372  INTEGER(KIND=large) :: ll
3373  INTEGER(KIND=large) :: lj
3374  INTEGER(KIND=large) :: indij
3375  INTEGER(KIND=large) :: indid
3376  INTEGER(KIND=large) :: ij
3377  INTEGER :: ichunk
3378  !$ INTEGER OMP_GET_THREAD_NUM
3379  SAVE
3380  ! ...
3381  !$ DO i=1,n
3382  !$ b(i)=0.0D0 ! reset 'global' B()
3383  !$ END DO
3384  ichunk=min((n+mthrd-1)/mthrd/8+1,1024)
3385  IF(matsto == 1) THEN
3386  ! full symmetric matrix
3387  ! parallelize row loop
3388  ! private copy of B(N) for each thread, combined at end, init with 0.
3389  ! slot of 1024 'I' for next idle thread
3390  !$OMP PARALLEL DO &
3391  !$OMP PRIVATE(J,IJ) &
3392  !$OMP REDUCTION(+:B) &
3393  !$OMP SCHEDULE(DYNAMIC,ichunk)
3394  DO i=1,n
3395  ij=i
3396  ij=(ij*ij-ij)/2
3397  b(i)=globalmatd(ij+i)*x(i)
3398  DO j=1,i-1
3399  b(j)=b(j)+globalmatd(ij+j)*x(i)
3400  b(i)=b(i)+globalmatd(ij+j)*x(j)
3401  END DO
3402  END DO
3403  !$OMP END PARALLEL DO
3404  ELSE
3405  ! sparse, compressed matrix
3406  IF(sparsematrixoffsets(2,1) /= n+1) stop 'AVPROD: mismatched vector and matrix'
3407  iencdb=nencdb
3408  iencdm=ishft(1,iencdb)-1
3409  ! parallelize row loop
3410  ! slot of 1024 'I' for next idle thread
3411  !$OMP PARALLEL DO &
3412  !$OMP PRIVATE(IR,K,KK,LL,KL,KU,INDID,INDIJ,J,JC,JN,LJ,JJ) &
3413  !$OMP REDUCTION(+:B) &
3414  !$OMP SCHEDULE(DYNAMIC,ichunk)
3415  DO i=1,n
3416  iproc=0
3417  !$ IPROC=OMP_GET_THREAD_NUM() ! thread number
3418  b(i)=globalmatd(i)*x(i) ! diagonal elements
3419  ! ! off-diagonals double precision
3420  ir=i
3421  kk=sparsematrixoffsets(1,ir) ! offset in 'd' (column lists)
3422  ll=sparsematrixoffsets(2,ir) ! offset in 'j' (matrix)
3423  kl=0
3424  ku=sparsematrixoffsets(1,ir+1)-1-kk
3425  indid=kk
3426  indij=ll
3427  IF (sparsematrixcolumns(indid) /= 0) THEN ! no compression
3428  DO k=kl,ku
3429  j=sparsematrixcolumns(indid+k)
3430  b(j)=b(j)+globalmatd(indij+k)*x(i)
3431  b(i)=b(i)+globalmatd(indij+k)*x(j)
3432  END DO
3433  ELSE
3434  lj=0
3435  ku=((ku+1)*8)/9-1 ! number of regions (-1)
3436  indid=indid+ku/8+1 ! skip group offsets
3437  DO kl=0,ku
3438  jc=sparsematrixcolumns(indid+kl)
3439  j=ishft(jc,-iencdb)
3440  jn=iand(jc, iencdm)
3441  DO jj=1,jn
3442  b(j)=b(j)+globalmatd(indij+lj)*x(i)
3443  b(i)=b(i)+globalmatd(indij+lj)*x(j)
3444  j=j+1
3445  lj=lj+1
3446  END DO
3447  END DO
3448  END IF
3449 
3450  IF (nspc > 1) THEN
3451  ir=i+n+1 ! off-diagonals single precision
3452  kk=sparsematrixoffsets(1,ir) ! offset in 'd' (column lists)
3453  ll=sparsematrixoffsets(2,ir) ! offset in '.' (matrix)
3454  kl=0
3455  ku=sparsematrixoffsets(1,ir+1)-1-kk
3456  indid=kk
3457  indij=ll
3458  IF (sparsematrixcolumns(indid) /= 0) THEN ! no compression
3459  DO k=kl,ku
3460  j=sparsematrixcolumns(indid+k)
3461  b(j)=b(j)+dble(globalmatf(indij+k))*x(i)
3462  b(i)=b(i)+dble(globalmatf(indij+k))*x(j)
3463  END DO
3464  ELSE
3465  lj=0
3466  ku=((ku+1)*8)/9-1 ! number of regions (-1)
3467  indid=indid+ku/8+1 ! skip group offsets
3468  DO kl=0,ku
3469  jc=sparsematrixcolumns(indid+kl)
3470  j=ishft(jc,-iencdb)
3471  jn=iand(jc, iencdm)
3472  DO jj=1,jn
3473  b(j)=b(j)+dble(globalmatf(indij+lj))*x(i)
3474  b(i)=b(i)+dble(globalmatf(indij+lj))*x(j)
3475  j=j+1
3476  lj=lj+1
3477  END DO
3478  END DO
3479  END IF
3480  END IF
3481  END DO
3482  !$OMP END PARALLEL DO
3483  ENDIF
3484 
3485 END SUBROUTINE avprod
3486 
3494 
3495 FUNCTION ijadd(itema,itemb) ! index using "d" and "z"
3496  USE mpmod
3497 
3498  IMPLICIT NONE
3499  INTEGER :: iencdb
3500  INTEGER :: iencdm
3501  INTEGER :: isgn
3502  INTEGER :: ispc
3503  INTEGER :: item2
3504  INTEGER :: jtem
3505  INTEGER :: jtemc
3506  INTEGER :: jtemn
3507 
3508  INTEGER, INTENT(IN) :: itema
3509  INTEGER, INTENT(IN) :: itemb
3510 
3511  INTEGER(KIND=large) :: ijadd
3512  INTEGER(KIND=large) :: k
3513  INTEGER(KIND=large) :: kk
3514  INTEGER(KIND=large) :: kl
3515  INTEGER(KIND=large) :: ku
3516  INTEGER(KIND=large) :: indid
3517  INTEGER(KIND=large) :: nd
3518  INTEGER(KIND=large) :: ll
3519  INTEGER(KIND=large) :: k8
3520  INTEGER(KIND=large) :: item1
3521  ! ...
3522  ijadd=0
3523  nd=sparsematrixoffsets(2,1)-1 ! dimension of matrix
3524  item1=max(itema,itemb) ! larger index
3525  item2=min(itema,itemb) ! smaller index
3526  IF(item2 <= 0.OR.item1 > nd) return
3527  IF(item1 == item2) THEN ! diagonal element
3528  ijadd=item1
3529  return
3530  END IF
3531  ! ! off-diagonal element
3532  iencdb=nencdb ! encoding info
3533  iencdm=ishft(1,iencdb)-1
3534  isgn=-1
3535  outer: DO ispc=1,nspc
3536  kk=sparsematrixoffsets(1,item1) ! offset in 'd' (column lists)
3537  ll=sparsematrixoffsets(2,item1) ! offset in 'j' (matrix)
3538  kl=0
3539  ku=sparsematrixoffsets(1,item1+1)-1-kk
3540  indid=kk
3541  item1=item1+nd+1
3542  isgn=-isgn
3543  IF (sparsematrixcolumns(indid) == 0) THEN ! compression ?
3544 
3545  ku=((ku+1)*8)/9-1 ! number of regions (-1)
3546  indid=indid+ku/8+1 ! skip group offsets
3547  kl=0
3548  IF(ku < kl) cycle outer ! not found
3549  DO
3550  k=(kl+ku)/2 ! binary search
3551  jtemc=sparsematrixcolumns(indid+k) ! compressed information
3552  jtem =ishft(jtemc,-iencdb) ! first column of region
3553  jtemn=jtem+iand(jtemc,iencdm) ! first column after region
3554  IF(item2 >= jtem.AND.item2 < jtemn) exit ! found
3555  IF(item2 < jtem) THEN
3556  ku=k-1
3557  ELSE IF(item2 >= jtemn) THEN
3558  kl=k+1
3559  END IF
3560  IF(kl <= ku) cycle
3561  cycle outer ! not found
3562  END DO
3563  k8=k/8 ! region group (-1)
3564  ll=ll+sparsematrixcolumns(kk+k8) ! offset for group of (8) regions
3565  DO kl=k8*8,k-1
3566  ll=ll+iand(sparsematrixcolumns(indid+kl),iencdm) ! add region lengths
3567  END DO
3568  ijadd=ll+item2-jtem
3569 
3570  ELSE
3571 
3572  IF(ku < kl) cycle outer ! not found
3573  DO
3574  k=(kl+ku)/2 ! binary search
3575  jtem=sparsematrixcolumns(indid+k)
3576  jtemn=jtem
3577  IF(item2 == jtem) exit ! found
3578  IF(item2 < jtem) THEN
3579  ku=k-1
3580  ELSE IF(item2 > jtem) THEN
3581  kl=k+1
3582  END IF
3583  IF(kl <= ku) cycle
3584  cycle outer ! not found
3585  END DO
3586  ijadd=ll+k
3587 
3588  END IF
3589  ijadd=ijadd*isgn
3590  return
3591  END DO outer
3592 
3593 END FUNCTION ijadd
3594 
3603 
3604 SUBROUTINE sechms(deltat,nhour,minut,secnd)
3605  IMPLICIT NONE
3606  REAL, INTENT(IN) :: deltat
3607  INTEGER, INTENT(OUT) :: minut
3608  INTEGER, INTENT(OUT):: nhour
3609  REAL, INTENT(OUT):: secnd
3610  INTEGER:: nsecd
3611  ! DELTAT = time in sec -> NHOUR,MINUT,SECND
3612  ! ...
3613  nsecd=nint(deltat) ! -> integer
3614  nhour=nsecd/3600
3615  minut=nsecd/60-60*nhour
3616  secnd=deltat-60*(minut+60*nhour)
3617 END SUBROUTINE sechms
3618 
3646 
3647 INTEGER FUNCTION inone(item) ! translate 1-D identifier to nrs
3648  USE mpmod
3649  USE mpdalc
3650 
3651  IMPLICIT NONE
3652  INTEGER, INTENT(IN) :: item
3653  INTEGER :: j
3654  INTEGER :: k
3655  INTEGER :: iprime
3656  INTEGER(KIND=large) :: length
3657  INTEGER(KIND=large), PARAMETER :: two = 2
3658 
3659  inone=0
3660  IF(item <= 0) return
3661  IF(globalparheader(-1) == 0) THEN
3662  length=128 ! initial number
3663  CALL mpalloc(globalparlabelindex,two,length,'INONE: label & index')
3664  CALL mpalloc(globalparhashtable,2*length,'INONE: hash pointer')
3665  globalparhashtable = 0
3666  globalparheader(-0)=int(length) ! length of labels/indices
3667  globalparheader(-1)=0 ! number of stored items
3668  globalparheader(-2)=0 ! =0 during build-up
3669  globalparheader(-3)=int(length) ! next number
3670  globalparheader(-4)=iprime(globalparheader(-0)) ! prime number
3671  globalparheader(-5)=0 ! number of overflows
3672  globalparheader(-6)=0 ! nr of variable parameters
3673  END IF
3674  outer: DO
3675  j=1+mod(item,globalparheader(-4))+globalparheader(-0)
3676  inner: DO ! normal case: find item
3677  k=j
3678  j=globalparhashtable(k)
3679  IF(j == 0) exit inner ! unused hash code
3680  IF(item == globalparlabelindex(1,j)) exit outer ! found
3681  END DO inner
3682  ! not found
3683  IF(globalparheader(-1) == globalparheader(-0).OR.globalparheader(-2) /= 0) THEN
3684  globalparheader(-5)=globalparheader(-5)+1 ! overflow
3685  j=0
3686  return
3687  END IF
3688  globalparheader(-1)=globalparheader(-1)+1 ! increase number of elements
3689  globalparheader(-3)=globalparheader(-1)
3690  j=globalparheader(-1)
3691  globalparhashtable(k)=j ! hash index
3692  globalparlabelindex(1,j)=item ! add new item
3693  globalparlabelindex(2,j)=0 ! reset counter
3694  IF(globalparheader(-1) /= globalparheader(-0)) exit outer
3695  ! update with larger dimension and redefine index
3696  globalparheader(-3)=globalparheader(-3)*2
3697  CALL upone
3698  IF (lvllog > 1) WRITE(lunlog,*) 'INONE: array increased to', &
3699  globalparheader(-3),' words'
3700  END DO outer
3701 
3702  IF(globalparheader(-2) == 0) THEN
3703  globalparlabelindex(2,j)=globalparlabelindex(2,j)+1 ! increase counter
3704  globalparheader(-7)=globalparheader(-7)+1
3705  END IF
3706  inone=j
3707 END FUNCTION inone
3708 
3710 SUBROUTINE upone
3711  USE mpmod
3712  USE mpdalc
3713 
3714  IMPLICIT NONE
3715  INTEGER :: i
3716  INTEGER :: j
3717  INTEGER :: k
3718  INTEGER :: iprime
3719  INTEGER :: nused
3720  LOGICAL :: finalupdate
3721  INTEGER(KIND=large) :: oldlength
3722  INTEGER(KIND=large) :: newlength
3723  INTEGER(KIND=large), PARAMETER :: two = 2
3724  INTEGER, DIMENSION(:,:), ALLOCATABLE :: temparr
3725  SAVE
3726  ! ...
3727  finalupdate=(globalparheader(-3) == globalparheader(-1))
3728  IF(finalupdate) THEN ! final (cleanup) call
3729  CALL sort2k(globalparlabelindex,globalparheader(-1)) ! sort items
3730  END IF
3731  ! save old LabelIndex
3732  nused = globalparheader(-1)
3733  oldlength = globalparheader(-0)
3734  CALL mpalloc(temparr,two,oldlength,'INONE: temp array')
3735  temparr(:,1:nused)=globalparlabelindex(:,1:nused)
3736  CALL mpdealloc(globalparlabelindex)
3737  CALL mpdealloc(globalparhashtable)
3738  ! create new LabelIndex
3739  newlength = globalparheader(-3)
3740  CALL mpalloc(globalparlabelindex,two,newlength,'INONE: label & index')
3741  CALL mpalloc(globalparhashtable,2*newlength,'INONE: hash pointer')
3742  globalparhashtable = 0
3743  globalparlabelindex(:,1:nused) = temparr(:,1:nused) ! copy back saved content
3744  CALL mpdealloc(temparr)
3745  globalparheader(-0)=int(newlength) ! length of labels/indices
3746  globalparheader(-3)=globalparheader(-1)
3747  globalparheader(-4)=iprime(globalparheader(-0)) ! prime number < LNDA
3748  ! redefine hash
3749  outer: DO i=1,globalparheader(-1)
3750  j=1+mod(globalparlabelindex(1,i),globalparheader(-4))+globalparheader(-0)
3751  inner: DO
3752  k=j
3753  j=globalparhashtable(k)
3754  IF(j == 0) exit inner ! unused hash code
3755  IF(j == i) cycle outer ! found
3756  ENDDO inner
3757  globalparhashtable(k)=i
3758  END DO outer
3759  IF(.NOT.finalupdate) return
3760 
3761  globalparheader(-2)=1 ! set flag to inhibit further updates
3762  IF (lvllog > 1) THEN
3763  WRITE(lunlog,*) ' '
3764  WRITE(lunlog,*) 'INONE: array reduced to',newlength,' words'
3765  WRITE(lunlog,*) 'INONE:',globalparheader(-1),' items stored.'
3766  END IF
3767 END SUBROUTINE upone ! update, redefine
3768 
3773 
3774 INTEGER FUNCTION iprime(n)
3775  IMPLICIT NONE
3776  INTEGER, INTENT(IN) :: n
3777  INTEGER :: nprime
3778  INTEGER :: nsqrt
3779  INTEGER :: i
3780  ! ...
3781  SAVE
3782  nprime=n ! max number
3783  IF(mod(nprime,2) == 0) nprime=nprime+1 ! ... odd number
3784  outer: DO
3785  nprime=nprime-2 ! next lower odd number
3786  nsqrt=ifix(sqrt(float(nprime)))
3787  DO i=3,nsqrt,2 !
3788  IF(i*(nprime/i) == nprime) cycle outer ! test prime number
3789  END DO
3790  exit outer ! found
3791  END DO outer
3792  iprime=nprime
3793 END FUNCTION iprime
3794 
3804 SUBROUTINE loop1
3805  USE mpmod
3806  USE mpdalc
3807 
3808  IMPLICIT NONE
3809  INTEGER :: i
3810  INTEGER :: idum
3811  INTEGER :: in
3812  INTEGER :: indab
3813  INTEGER :: itgbi
3814  INTEGER :: itgbl
3815  INTEGER :: ivgbi
3816  INTEGER :: j
3817  INTEGER :: mqi
3818  INTEGER :: nc21
3819  INTEGER :: nr
3820  INTEGER :: nwrd
3821  INTEGER :: inone
3822  REAL :: param
3823  REAL :: presg
3824  REAL :: prewt
3825 
3826  REAL :: plvs(3) ! vector array: real and ...
3827  INTEGER :: lpvs(3) ! ... integer
3828  equivalence(plvs(1),lpvs(1))
3829  INTEGER (KIND=large) :: length
3830  SAVE
3831  ! ...
3832  WRITE(lunlog,*) ' '
3833  WRITE(lunlog,*) 'LOOP1: starting'
3834  CALL mstart('LOOP1')
3835  ! add labels from parameter, constraints, measurements -------------
3836  DO i=1, lenparameters
3837  idum=inone(listparameters(i)%label)
3838  END DO
3839  DO i=1, lenpresigmas
3840  idum=inone(listpresigmas(i)%label)
3841  END DO
3842  DO i=1, lenconstraints
3843  idum=inone(listconstraints(i)%label)
3844  END DO
3845  DO i=1, lenmeasurements
3846  idum=inone(listmeasurements(i)%label)
3847  END DO
3848 
3849  IF(globalparheader(-1) /= 0) THEN
3850  WRITE(lunlog,*) 'LOOP1:',globalparheader(-1), ' labels from txt data stored'
3851  END IF
3852  WRITE(lunlog,*) 'LOOP1: reading data files'
3853 
3854  DO
3855  DO j=1,globalparheader(-1)
3856  globalparlabelindex(2,j)=0 ! reset count
3857  END DO
3858 
3859  ! read all data files and add all labels to global labels table ----
3860 
3861  IF(mprint /= 0) THEN
3862  WRITE(*,*) 'Read all binary data files:'
3863  END IF
3864  CALL hmpldf(1,'Number of words/record in binary file')
3865  CALL hmpdef(8,0.0,60.0,'not_stored data per record')
3866  ! define read buffer
3867  nc21=ncache/(21*mthrdr) ! split read cache 1 : 10 : 10 for pointers, ints, floats
3868  nwrd=nc21+1
3869  length=nwrd*mthrdr
3870  CALL mpalloc(readbufferpointer,length,'read buffer, pointer')
3871  nwrd=nc21*10+2+ndimbuf
3872  length=nwrd*mthrdr
3873  CALL mpalloc(readbufferdatai,length,'read buffer, integer')
3874  CALL mpalloc(readbufferdataf,length,'read buffer, float')
3875 
3876  DO
3877  CALL peread(nr) ! read records
3878  IF (skippedrecords == 0) CALL peprep(0) ! prepare records
3879  IF(nr <= 0) exit ! end of data?
3880  END DO
3881  ! release read buffer
3882  CALL mpdealloc(readbufferdataf)
3883  CALL mpdealloc(readbufferdatai)
3884  CALL mpdealloc(readbufferpointer)
3885  IF (skippedrecords == 0) THEN
3886  exit
3887  ELSE
3888  WRITE(lunlog,*) 'LOOP1: reading data files again'
3889  END IF
3890  END DO
3891 
3892  IF(nhistp /= 0) THEN
3893  CALL hmprnt(1)
3894  CALL hmprnt(8)
3895  END IF
3896  CALL hmpwrt(1)
3897  CALL hmpwrt(8)
3898  CALL upone ! finalize the global label table
3899  ntgb = globalparheader(-1) ! total number of labels/parameters
3900  WRITE(lunlog,*) 'LOOP1:',ntgb, &
3901  ' is total number NTGB of labels/parameters'
3902  ! histogram number of entries per label ----------------------------
3903  CALL hmpldf(2,'Number of entries per label')
3904  DO j=1,ntgb
3905  CALL hmplnt(2,globalparlabelindex(2,j))
3906  END DO
3907  IF(nhistp /= 0) CALL hmprnt(2) ! print histogram
3908  CALL hmpwrt(2) ! write to his file
3909 
3910  ! three subarrays for all global parameters ------------------------
3911  length=ntgb
3912  CALL mpalloc(globalparameter,length,'global parameters')
3913  globalparameter=0.0d0
3914  CALL mpalloc(globalparpresigma,length,'pre-sigmas') ! presigmas
3915  globalparpresigma=0.
3916  CALL mpalloc(globalparstart,length,'global parameters at start')
3917  globalparstart=0.
3918  CALL mpalloc(globalparcopy,length,'copy of global parameters')
3919 
3920  DO i=1,lenparameters ! parameter start values
3921  param=listparameters(i)%value
3922  in=inone(listparameters(i)%label)
3923  IF(in /= 0) THEN
3924  globalparameter(in)=param
3925  globalparstart(in)=param
3926  ENDIF
3927  END DO
3928 
3929  npresg=0
3930  DO i=1,lenpresigmas ! pre-sigma values
3931  presg=listpresigmas(i)%value
3932  in=inone(listpresigmas(i)%label)
3933  IF(in /= 0) THEN
3934  IF(presg > 0.0) npresg=npresg+1 ! FIXME: check if enough 'entries'?
3935  globalparpresigma(in)=presg ! insert pre-sigma 0 or > 0
3936  END IF
3937  END DO
3938  WRITE(lunlog,*) 'LOOP1:',npresg,' is number of pre-sigmas'
3939  WRITE(*,*) 'LOOP1:',npresg,' is number of pre-sigmas'
3940  IF(npresg == 0) WRITE(*,*) 'Warning: no pre-sigmas defined'
3941 
3942  ! determine flag variable (active) or fixed (inactive) -------------
3943 
3944  indab=0
3945  DO i=1,ntgb
3946  IF(globalparlabelindex(2,i) >= mreqen.AND.globalparpresigma(i) >= 0.0) THEN
3947  indab=indab+1
3948  globalparlabelindex(2,i)=indab ! variable, used in matrix (active)
3949  ELSE
3950  globalparlabelindex(2,i)=-1 ! fixed, not used in matrix (not active)
3951  END IF
3952  END DO
3953  globalparheader(-6)=indab ! counted variable
3954  nvgb=indab ! nr of variable parameters
3955  WRITE(lunlog,*) 'LOOP1:',nvgb, ' is number NVGB of variable parameters'
3956 
3957  ! translation table of length NVGB of total global indices ---------
3958  length=nvgb
3959  CALL mpalloc(globalparvartototal,length,'translation table var -> total')
3960  indab=0
3961  DO i=1,ntgb
3962  IF(globalparlabelindex(2,i) > 0) THEN
3963  indab=indab+1
3964  globalparvartototal(indab)=i
3965  END IF
3966  END DO
3967 
3968  ! regularization ---------------------------------------------------
3969  CALL mpalloc(globalparpreweight,length,'pre-sigmas weights') ! presigma weights
3970  WRITE(*,112) ' Default pre-sigma =',regpre, &
3971  ' (if no individual pre-sigma defined)'
3972  WRITE(*,*) 'Pre-sigma factor is',regula
3973 
3974  IF(nregul == 0) THEN
3975  WRITE(*,*) 'No regularization will be done'
3976  ELSE
3977  WRITE(*,*) 'Regularization will be done, using factor',regula
3978  END IF
3979 112 format(a,e9.2,a)
3980  IF (nvgb <= 0) stop '... no variable global parameters'
3981 
3982  DO ivgbi=1,nvgb ! IVGBI = variable parameter index
3983  itgbi=globalparvartototal(ivgbi) ! ITGBI = global parameter index
3984  presg=globalparpresigma(itgbi) ! get pre-sigma
3985  prewt=0.0 ! pre-weight
3986  IF(presg > 0.0) THEN
3987  prewt=1.0/presg**2 ! 1/presigma^2
3988  ELSE IF(presg == 0.0.AND.regpre > 0.0) THEN
3989  prewt=1.0/regpre**2 ! default 1/presigma^2
3990  END IF
3991  globalparpreweight(ivgbi)=regula*prewt ! weight = factor / presigma^2
3992  END DO
3993 
3994  ! WRITE(*,*) 'GlPa_index GlPa_label array1 array6'
3995  DO i=1,ntgb
3996  itgbl=globalparlabelindex(1,i)
3997  ivgbi=globalparlabelindex(2,i)
3998  IF(ivgbi > 0) THEN
3999  ! WRITE(*,111) I,ITGBL,QM(IND1+I),QM(IND6+IVGBI)
4000  ELSE
4001  ! WRITE(*,111) I,ITGBL,QM(IND1+I)
4002  END IF
4003  END DO
4004  ! 111 FORMAT(I5,I10,F10.5,E12.4)
4005  WRITE(*,101) 'NTGB',ntgb,'total number of parameters'
4006  WRITE(*,101) 'NVGB',nvgb,'number of variable parameters'
4007 
4008  ! print overview over important numbers ----------------------------
4009 
4010  nrecal=nrec
4011  IF(mprint /= 0) THEN
4012  WRITE(*,*) ' '
4013  WRITE(*,101) ' NREC',nrec,'number of records'
4014  WRITE(*,101) 'MREQEN',mreqen,'required number of entries'
4015  IF (mreqpe > 1) WRITE(*,101) &
4016  'MREQPE',mreqpe,'required number of pair entries'
4017  IF (msngpe >= 1) WRITE(*,101) &
4018  'MSNGPE',msngpe,'max pair entries single prec. storage'
4019  WRITE(*,101) 'NTGB',ntgb,'total number of parameters'
4020  WRITE(*,101) 'NVGB',nvgb,'number of variable parameters'
4021  IF(mprint > 1) THEN
4022  WRITE(*,*) ' '
4023  WRITE(*,*) 'Global parameter labels:'
4024  mqi=ntgb
4025  IF(mqi <= 100) THEN
4026  WRITE(*,*) (globalparlabelindex(2,i),i=1,mqi)
4027  ELSE
4028  WRITE(*,*) (globalparlabelindex(2,i),i=1,30)
4029  WRITE(*,*) ' ...'
4030  mqi=((mqi-20)/20)*20+1
4031  WRITE(*,*) (globalparlabelindex(2,i),i=mqi,ntgb)
4032  END IF
4033  END IF
4034  WRITE(*,*) ' '
4035  WRITE(*,*) ' '
4036  END IF
4037  WRITE(8,*) ' '
4038  WRITE(8,101) ' NREC',nrec,'number of records'
4039  WRITE(8,101) 'MREQEN',mreqen,'required number of entries'
4040 
4041  WRITE(lunlog,*) 'LOOP1: ending'
4042  WRITE(lunlog,*) ' '
4043  CALL mend
4044 
4045 101 format(1x,a6,' =',i10,' = ',a)
4046 END SUBROUTINE loop1
4047 
4058 
4059 SUBROUTINE loop2
4060  USE mpmod
4061  USE mpdalc
4062 
4063  IMPLICIT NONE
4064  REAL :: chin2
4065  REAL :: chin3
4066  REAL :: cpr
4067  REAL :: fsum
4068  REAL :: gbc
4069  REAL :: gbu
4070  REAL :: glder
4071  INTEGER :: i
4072  INTEGER :: ia
4073  INTEGER :: ib
4074  INTEGER :: ibuf
4075  INTEGER :: icgb
4076  INTEGER :: iext
4077  INTEGER :: ihis
4078  INTEGER :: ij
4079  INTEGER :: ijn
4080  INTEGER :: inder
4081  INTEGER :: ioff
4082  INTEGER :: iproc
4083  INTEGER :: irecmm
4084  INTEGER :: isfrst
4085  INTEGER :: islast
4086  INTEGER :: ist
4087  INTEGER :: itgbi
4088  INTEGER :: itgbij
4089  INTEGER :: itgbik
4090  INTEGER :: ivgbij
4091  INTEGER :: ivgbik
4092  INTEGER :: j
4093  INTEGER :: ja
4094  INTEGER :: jb
4095  INTEGER :: jcmprs
4096  INTEGER :: jext
4097  INTEGER :: jsp
4098  INTEGER :: k
4099  INTEGER :: kfile
4100  INTEGER :: l
4101  INTEGER :: label
4102  INTEGER :: last
4103  INTEGER :: lu
4104  INTEGER :: lun
4105  INTEGER :: maeqnf
4106  INTEGER :: naeqna
4107  INTEGER :: naeqnf
4108  INTEGER :: naeqng
4109  INTEGER :: nc21
4110  INTEGER :: ncachd
4111  INTEGER :: ncachi
4112  INTEGER :: ncachr
4113  INTEGER :: nda
4114  INTEGER :: ndf
4115  INTEGER :: ndfmax
4116  INTEGER :: nfixed
4117  INTEGER :: nggd
4118  INTEGER :: nggi
4119  INTEGER :: nmatmo
4120  INTEGER :: noff
4121  INTEGER :: nr
4122  INTEGER :: nrecf
4123  INTEGER :: nrecmm
4124  INTEGER :: nst
4125  INTEGER :: nwrd
4126  INTEGER :: inone
4127  REAL :: wgh
4128  REAL :: wolfc3
4129  REAL :: wrec
4130  REAL :: chindl
4131 
4132  DOUBLE PRECISION::dstat(3)
4133  INTEGER*8:: noff8
4134  INTEGER(KIND=large):: ndimbi
4135  INTEGER(KIND=large):: ndimsa(4)
4136  INTEGER(KIND=large):: ndgn
4137  INTEGER(KIND=large):: matsiz(2)
4138  INTEGER(KIND=large):: matwords
4139  INTEGER(KIND=large):: length
4140  INTEGER(KIND=large):: rows
4141  INTEGER(KIND=large):: cols
4142  INTEGER(KIND=large), PARAMETER :: two=2
4143  INTEGER:: maxglobalpar = 0
4144  INTEGER:: maxlocalpar = 0
4145  INTEGER:: maxequations = 0
4146 
4147  INTERFACE ! needed for assumed-shape dummy arguments
4148  SUBROUTINE ndbits(ndims,ncmprs,nsparr,mnpair,ihst,jcmprs)
4149  USE mpdef
4150  INTEGER(KIND=large), DIMENSION(4), INTENT(OUT) :: ndims
4151  INTEGER, DIMENSION(:), INTENT(OUT) :: ncmprs
4152  INTEGER(KIND=large), DIMENSION(:,:), INTENT(OUT) :: nsparr
4153  INTEGER, INTENT(IN) :: mnpair
4154  INTEGER, INTENT(IN) :: ihst
4155  INTEGER, INTENT(IN) :: jcmprs
4156  END SUBROUTINE ndbits
4157  SUBROUTINE spbits(nsparr,nsparc,ncmprs)
4158  USE mpdef
4159  INTEGER(KIND=large), DIMENSION(:,:), INTENT(IN) :: nsparr
4160  INTEGER, DIMENSION(:), INTENT(OUT) :: nsparc
4161  INTEGER, DIMENSION(:), INTENT(IN) :: ncmprs
4162  END SUBROUTINE spbits
4163  END INTERFACE
4164 
4165  SAVE
4166 
4167  !$ INTEGER :: OMP_GET_THREAD_NUM
4168 
4169  isfrst(ibuf)=readbufferpointer(ibuf)+1
4170  islast(ibuf)=readbufferdatai(readbufferpointer(ibuf))
4171  inder(i)=readbufferdatai(i)
4172  glder(i)=readbufferdataf(i)
4173  ! ...
4174  WRITE(lunlog,*) ' '
4175  WRITE(lunlog,*) 'LOOP2: starting'
4176  CALL mstart('LOOP2')
4177 
4178  ! two subarrays to get the global parameter indices, used in an event
4179  length=nvgb
4180  CALL mpalloc(globalindexusage,length,'global index')
4181  CALL mpalloc(backindexusage,length,'back index')
4182  backindexusage=0
4183 
4184  ! constraints - determine number of constraints NCGB
4185  ncgb=0
4186  i=0
4187  last=-1
4188  ! find next constraint header and count nr of constraints
4189  DO WHILE(i < lenconstraints)
4190  i=i+1
4191  label=listconstraints(i)%label
4192  IF(last == 0.AND.label == (-1)) ncgb=ncgb+1
4193  last=label
4194  END DO
4195  WRITE(*,*) 'LOOP2:',ncgb,' constraints'
4196 
4197  nagb=nvgb+ncgb ! total number of fit parameters
4198  noff8=int8(nagb)*int8(nagb-1)/2
4199 
4200  ! read all data files and add all variable index pairs -------------
4201 
4202  IF(matsto == 2) THEN
4203  IF (mcmprs /= 0) numbit=max(numbit,2) ! identify single entries for compression
4204  CALL clbits(nagb,ndimbi,nencdb,numbit) ! get dimension for bit storage, encoding info
4205  END IF
4206 
4207  ! reading events===reading events===reading events===reading events=
4208  nrecf =0 ! records with fixed global parameters
4209  naeqng=0 ! count number of equations (with global der.)
4210  naeqnf=0 ! count number of equations ( " , fixed)
4211  naeqna=0 ! all
4212  WRITE(lunlog,*) 'LOOP2: start event reading'
4213  ! monitoring for sparse matrix?
4214  irecmm=0
4215  IF (matsto == 2.AND.matmon /= 0) THEN
4216  nmatmo=0
4217  IF (matmon > 0) THEN
4218  nrecmm=matmon
4219  ELSE
4220  nrecmm=1
4221  END IF
4222  END IF
4223  DO k=1,3
4224  dstat(k)=0.0d0
4225  END DO
4226  ! define read buffer
4227  nc21=ncache/(21*mthrdr) ! split read cache 1 : 10 : 10 for pointers, ints, floats
4228  nwrd=nc21+1
4229  length=nwrd*mthrdr
4230  CALL mpalloc(readbufferpointer,length,'read buffer, pointer')
4231  nwrd=nc21*10+2+ndimbuf
4232  length=nwrd*mthrdr
4233  CALL mpalloc(readbufferdatai,length,'read buffer, integer')
4234  CALL mpalloc(readbufferdataf,length,'read buffer, float')
4235 
4236  DO
4237  CALL peread(nr) ! read records
4238  CALL peprep(1) ! prepare records
4239  ioff=0
4240  DO ibuf=1,numreadbuffer ! buffer for current record
4241  nrec=readbufferdatai(isfrst(ibuf)-2) ! record
4242  ! Printout for DEBUG
4243  IF(nrec <= mdebug) THEN
4244  nda=0
4245  kfile=nint(readbufferdataf(isfrst(ibuf)-1)) ! file
4246  wrec =readbufferdataf(isfrst(ibuf)-2) ! weight
4247  WRITE(*,*) ' '
4248  WRITE(*,*) 'Record number ',nrec,' from file ',kfile
4249  IF (wgh /= 1.0) WRITE(*,*) ' weight ',wrec
4250  ist=isfrst(ibuf)
4251  nst=islast(ibuf)
4252  DO ! loop over measurements
4253  CALL isjajb(nst,ist,ja,jb,jsp)
4254  IF(ja == 0) exit
4255  nda=nda+1
4256  IF(nda > mdebg2) THEN
4257  IF(nda == mdebg2+1) WRITE(*,*) '... and more data'
4258  cycle
4259  END IF
4260  WRITE(*,*) ' '
4261  WRITE(*,*) nda, ' Measured value =',glder(ja),' +- ',glder(jb)
4262  WRITE(*,*) 'Local derivatives:'
4263  WRITE(*,107) (inder(ja+j),glder(ja+j),j=1,jb-ja-1)
4264 107 format(6(i3,g12.4))
4265  IF (jb < ist) THEN
4266  WRITE(*,*) 'Global derivatives:'
4267  WRITE(*,108) (globalparlabelindex(1,inder(jb+j)),inder(jb+j), &
4268  globalparlabelindex(2,inder(jb+j)),glder(jb+j),j=1,ist-jb)
4269 108 format(3i11,g12.4)
4270  END IF
4271  IF(nda == 1) THEN
4272  WRITE(*,*) 'total_par_label __label__ var_par_index derivative'
4273  END IF
4274  END DO
4275  WRITE(*,*) ' '
4276  END IF
4277 
4278  nagbn =0 ! count number of global derivatives
4279  nalcn =0 ! count number of local derivatives
4280  naeqn =0 ! count number of equations
4281  maeqnf=naeqnf
4282  ist=isfrst(ibuf)
4283  nst=islast(ibuf)
4284  nwrd=nst-ist+1
4285  DO ! loop over measurements
4286  CALL isjajb(nst,ist,ja,jb,jsp)
4287  IF(ja == 0.AND.jb == 0) exit
4288  naeqn=naeqn+1
4289  naeqna=naeqna+1
4290  IF(ja /= 0) THEN
4291  IF (ist > jb) naeqng=naeqng+1
4292  nfixed=0
4293  DO j=1,ist-jb
4294  ij=inder(jb+j) ! index of global parameter
4295  ij=globalparlabelindex(2,ij) ! change to variable parameter
4296  IF(ij > 0) THEN
4297  ijn=backindexusage(ij) ! get index of index
4298  IF(ijn == 0) THEN ! not yet included
4299  nagbn=nagbn+1 ! count
4300  globalindexusage(nagbn)=ij ! store variable index
4301  backindexusage(ij)=nagbn ! store back index
4302  END IF
4303  ELSE
4304  nfixed=nfixed+1
4305  END IF
4306  END DO
4307  IF (nfixed > 0) naeqnf=naeqnf+1
4308  END IF
4309 
4310  IF(ja /= 0.AND.jb /= 0) THEN
4311  DO j=1,jb-ja-1 ! local parameters
4312  ij=inder(ja+j)
4313  nalcn=max(nalcn,ij)
4314  END DO
4315  END IF
4316  END DO
4317 
4318  ! end-of-event
4319  IF (naeqnf > maeqnf) nrecf=nrecf+1
4320  irecmm=irecmm+1
4321  ! end-of-event-end-of-event-end-of-event-end-of-event-end-of-event-e
4322 
4323  maxglobalpar=max(nagbn,maxglobalpar) ! maximum number of global parameters
4324  maxlocalpar=max(nalcn,maxlocalpar) ! maximum number of local parameters
4325  maxequations=max(naeqn,maxequations) ! maximum number of equations
4326 
4327  ! sample statistics for caching
4328  dstat(1)=dstat(1)+dfloat((nwrd+2)*2) ! record size
4329  dstat(2)=dstat(2)+dfloat(nagbn+2) ! indices,
4330  dstat(3)=dstat(3)+dfloat(nagbn*nagbn+nagbn) ! data for MUPDAT
4331 
4332  CALL sort1k(globalindexusage,nagbn) ! sort global par.
4333  ! overwrite read buffer with lists of global labels
4334  ioff=ioff+1
4335  readbufferpointer(ibuf)=ioff
4336  readbufferdatai(ioff)=ioff+nagbn
4337  DO i=1,nagbn ! reset global index array
4338  iext=globalindexusage(i)
4339  backindexusage(iext)=0
4340  readbufferdatai(ioff+i)=iext
4341  END DO
4342  ioff=ioff+nagbn
4343 
4344  END DO
4345  ioff=0
4346 
4347  IF (matsto == 2) THEN
4348  !$OMP PARALLEL &
4349  !$OMP DEFAULT(PRIVATE) &
4350  !$OMP SHARED(numReadBuffer,readBufferPointer,readBufferDataI,MTHRD)
4351  iproc=0
4352  !$ IPROC=OMP_GET_THREAD_NUM() ! thread number
4353  DO ibuf=1,numreadbuffer
4354  ist=isfrst(ibuf)
4355  nst=islast(ibuf)
4356  DO i=ist,nst ! store all combinations
4357  iext=readbufferdatai(i) ! variable global index
4358  !$ IF (MOD(IEXT,MTHRD).EQ.IPROC) THEN ! distinct rows per thread
4359  DO l=ist,i
4360  jext=readbufferdatai(l)
4361  CALL inbits(iext,jext,1) ! save space
4362  END DO
4363  !$ ENDIF
4364  END DO
4365  END DO
4366  !$OMP END PARALLEL
4367  ! monitoring
4368  IF (matmon /= 0.AND. &
4369  (irecmm >= nrecmm.OR.irecmm == mxrec)) THEN
4370  IF (nmatmo == 0) THEN
4371  WRITE(*,*)
4372  WRITE(*,*) 'Monitoring of sparse matrix construction'
4373  WRITE(*,*) ' records ........ off-diagonal elements ', &
4374  '....... compression memory'
4375  WRITE(*,*) ' non-zero used(double) used', &
4376  '(float) [%] [GB]'
4377  END IF
4378  nmatmo=nmatmo+1
4379  jcmprs=max(mcmprs,msngpe)
4380  CALL ckbits(ndimsa,mreqpe,jcmprs)
4381  gbc=4.0e-9*float(ndimsa(2)+2*ndimsa(3)+ndimsa(4)) ! GB compressed
4382  gbu=1.2e-8*float(ndimsa(3)+ndimsa(4)) ! GB uncompressed
4383  cpr=100.0*gbc/gbu
4384  WRITE(*,1177) irecmm,ndimsa(1),ndimsa(3),ndimsa(4),cpr,gbc
4385 1177 format(i9,3i13,f10.2,f11.6)
4386  DO WHILE(irecmm >= nrecmm)
4387  IF (matmon > 0) THEN
4388  nrecmm=nrecmm+matmon
4389  ELSE
4390  nrecmm=nrecmm*2
4391  END IF
4392  END DO
4393  END IF
4394 
4395  END IF
4396 
4397  IF (nr <= 0) exit ! next block of events ?
4398  END DO
4399  ! release read buffer
4400  CALL mpdealloc(readbufferdataf)
4401  CALL mpdealloc(readbufferdatai)
4402  CALL mpdealloc(readbufferpointer)
4403 
4404  WRITE(lunlog,*) 'LOOP2: event reading ended - end of data'
4405  DO k=1,3
4406  dstat(k)=dstat(k)/dfloat(nrec)
4407  END DO
4408  ! end=of=data=end=of=data=end=of=data=end=of=data=end=of=data=end=of
4409 
4410  IF(matsto == 2) THEN
4411 
4412  ! constraints and index pairs with Lagrange multiplier
4413 
4414 
4415  ! constraints - determine number of constraints NCGB and index-pairs
4416  ! Lagrange multiplier and global parameters
4417 
4418 
4419  i=0
4420  icgb=0
4421  last=-1
4422  ! find next constraint header
4423  DO WHILE(i < lenconstraints)
4424  i=i+1
4425  label=listconstraints(i)%label
4426  IF(last == 0.AND.label == (-1)) icgb=icgb+1
4427  IF(label > 0) THEN
4428  itgbi=inone(label)
4429  ij=globalparlabelindex(2,itgbi) ! change to variable parameter
4430  IF(ij > 0) THEN
4431  CALL inbits(nvgb+icgb,ij,mreqpe)
4432  END IF
4433  END IF
4434  last=label
4435  END DO
4436 
4437  ! measurements - determine index-pairs
4438 
4439 
4440  i=1
4441  DO WHILE (i <= lenmeasurements)
4442  i=i+2
4443  ! loop over label/factor pairs
4444  ia=i
4445  DO
4446  i=i+1
4447  IF(i > lenmeasurements) exit
4448  IF(listmeasurements(i)%label == 0) exit
4449  END DO
4450  ib=i-1
4451 
4452  DO j=ia,ib
4453  itgbij=inone(listmeasurements(j)%label) ! total parameter index
4454  ! first index
4455  ivgbij=0
4456  IF(itgbij /= 0) ivgbij=globalparlabelindex(2,itgbij) ! variable-parameter index
4457  DO k=ia,j
4458  itgbik=inone(listmeasurements(k)%label) ! total parameter index
4459  ! second index
4460  ivgbik=0
4461  IF(itgbik /= 0) ivgbik=globalparlabelindex(2,itgbik) ! variable-parameter index
4462  IF(ivgbij > 0.AND.ivgbik > 0) THEN
4463  CALL inbits(ivgbij,ivgbik,mreqpe)
4464  IF (mprint > 1) WRITE(*,*) 'add index pair ',ivgbij,ivgbik
4465  END IF
4466  END DO
4467  END DO
4468 
4469  END DO
4470  END IF
4471 
4472  ! print numbers ----------------------------------------------------
4473 
4474  IF (nagb >= 65536) THEN
4475  noff=int(noff8/1000)
4476  ELSE
4477  noff=int(noff8)
4478  END IF
4479  ndgn=0
4480  nspc=1 ! number of precision types (double, single) for matrix storage
4481  matwords=0
4482  IF(matsto == 2) THEN
4483  ihis=0
4484  IF (mhispe > 0) THEN
4485  ihis=15
4486  CALL hmpdef(ihis,0.0,float(mhispe), 'NDBITS: #off-diagonal elements')
4487  END IF
4488  jcmprs=max(mcmprs,msngpe)
4489  IF (jcmprs > 0.AND.numbit > 1) nspc=2 ! mixed precision storage
4490  length=nagb*nspc
4491  CALL mpalloc(sparsematrixcompression,length,'INBITS: row compression')
4492  sparsematrixcompression=0
4493  length=(nagb+1)*nspc
4494  CALL mpalloc(sparsematrixoffsets,two,length, 'sparse matrix row offsets')
4495  CALL ndbits(ndimsa,sparsematrixcompression,sparsematrixoffsets, &
4496  mreqpe,ihis,jcmprs)
4497  ndgn=ndimsa(3)+ndimsa(4) ! actual number of off-diagonal elements
4498  matwords=ndimsa(2)+length ! size of sparsity structure
4499 
4500  IF (mhispe > 0) THEN
4501  IF (nhistp /= 0) CALL hmprnt(ihis)
4502  CALL hmpwrt(ihis)
4503  END IF
4504  END IF
4505 
4506  nagbn=maxglobalpar ! max number of global parameters in one event
4507  nalcn=maxlocalpar ! max number of local parameters in one event
4508  naeqn=maxequations ! max number of equations in one event
4509  CALL mpdealloc(globalindexusage)
4510  CALL mpdealloc(backindexusage)
4511  ! matrices for event matrices
4512  ! split up cache
4513  IF (fcache(2) == 0.0) THEN ! from data (DSTAT)
4514  fcache(1)=sngl(dstat(1))*fcache(1) ! leave some part free for fluctuations
4515  fcache(2)=sngl(dstat(2))
4516  fcache(3)=sngl(dstat(3))
4517  END IF
4518  fsum=fcache(1)+fcache(2)+fcache(3)
4519  DO k=1,3
4520  fcache(k)=fcache(k)/fsum
4521  END DO
4522  ncachr=ifix(float(ncache)*fcache(1)+0.5) ! read cache
4523  ! define read buffer
4524  nc21=ncachr/(21*mthrdr) ! split read cache 1 : 10 : 10 for pointers, ints, floats
4525  nwrd=nc21+1
4526  length=nwrd*mthrdr
4527  CALL mpalloc(readbufferpointer,length,'read buffer, pointer')
4528  nwrd=nc21*10+2+ndimbuf
4529  length=nwrd*mthrdr
4530  CALL mpalloc(readbufferdatai,length,'read buffer, integer')
4531  CALL mpalloc(readbufferdataf,length,'read buffer, float')
4532 
4533  ncachi=ifix(float(ncache)*fcache(2)+0.5) ! index cache
4534  ncachd=ncache-ncachr-ncachi ! data cache
4535  nggd=(nagbn*nagbn+nagbn)/2+ncachd/(2*mthrd) ! number of double
4536  nggi=2+nagbn+ncachi/mthrd ! number of ints
4537  length=nagbn*mthrd
4538  CALL mpalloc(globalindexusage,length, 'global parameters (dim =max/event)')
4539  length=nvgb*mthrd
4540  CALL mpalloc(backindexusage,length,'global variable-index array')
4541  backindexusage=0
4542  length=nagbn*nalcn
4543  CALL mpalloc(localglobalmatrix,length,'local/global matrix')
4544  length=nggd*mthrd
4545  CALL mpalloc(writebufferupdates,length,'symmetric update matrices')
4546  writebufferheader(-1)=nggd ! number of words per thread
4547  writebufferheader(-2)=(nagbn*nagbn+nagbn)/2 ! min free (double) words
4548  length=nggi*mthrd
4549  CALL mpalloc(writebufferindices,length,'symmetric update matrix indices')
4550  rows=7; cols=mthrd
4551  CALL mpalloc(writebufferinfo,rows,cols,'write buffer status (I)')
4552  rows=2; cols=mthrd
4553  CALL mpalloc(writebufferdata,rows,cols,'write buffer status (F)')
4554  writebufferheader(1)=nggi ! number of words per thread
4555  writebufferheader(2)=nagbn+2 ! min free words
4556 
4557  ! print all relevant dimension parameters
4558 
4559  DO lu=6,8,2 ! unit 6 and 8
4560 
4561  WRITE(*,*) ' '
4562  WRITE(lu,101) 'NTGB',ntgb,'total number of parameters'
4563  WRITE(lu,102) '(all parameters, appearing in binary files)'
4564  WRITE(lu,101) 'NVGB',nvgb,'number of variable parameters'
4565  WRITE(lu,102) '(appearing in fit matrix/vectors)'
4566  WRITE(lu,101) 'NAGB',nagb,'number of fit parameters'
4567  WRITE(lu,102) '(including Lagrange multiplier or reduced)'
4568  WRITE(lu,101) 'MBANDW',mbandw,'band width of band matrix'
4569  WRITE(lu,102) '(if =0, no band matrix)'
4570  IF (nagb >= 65536) THEN
4571  WRITE(lu,101) 'NOFF/K',noff,'max number of off-diagonal elements'
4572  ELSE
4573  WRITE(lu,101) 'NOFF',noff,'max number of off-diagonal elements'
4574  END IF
4575  IF(ndgn /= 0) &
4576  WRITE(lu,101) 'NDGN',ndgn,'actual number of off-diagonal elements'
4577  WRITE(lu,101) 'NCGB',ncgb,'number of constraints'
4578  WRITE(lu,101) 'NAGBN',nagbn,'max number of global parameters in an event'
4579  WRITE(lu,101) 'NALCN',nalcn,'max number of local parameters in an event'
4580  WRITE(lu,101) 'NAEQN',naeqn,'max number of equations in an event'
4581  IF (mprint > 1) THEN
4582  WRITE(lu,101) 'NAEQNA',naeqna,'number of equations'
4583  WRITE(lu,101) 'NAEQNG',naeqng, &
4584  'number of equations with global derivatives'
4585  WRITE(lu,101) 'NAEQNF',naeqnf, &
4586  'number of equations with fixed global derivatives'
4587  WRITE(lu,101) 'NRECF',nrecf, &
4588  'number of records with fixed global derivatives'
4589  END IF
4590  IF (ncache > 0) THEN
4591  WRITE(lu,101) 'NCACHE',ncache,'number of words for caching'
4592  WRITE(lu,111) (fcache(k)*100.0,k=1,3)
4593 111 format(22x,'cache splitting ',3(f6.1,' %'))
4594  END IF
4595  WRITE(lu,*) ' '
4596  IF(matsto == 2) THEN
4597  ! FRUSE=FLOAT(MQ(INDF-1))/FLOAT(NOFF)
4598  ! WRITE(LU,123) 100.0*FRUSE
4599  END IF
4600  ! 123 FORMAT(' Fraction of non-zero off-diagonal elements is',
4601  ! + F6.1,' %')
4602 
4603  ! print information about methods and matrix storage modes
4604 
4605  WRITE(lu,*) ' '
4606  WRITE(lu,*) 'Solution method and matrix-storage mode:'
4607  IF(metsol == 1) THEN
4608  WRITE(lu,*) ' METSOL = 1: matrix inversion'
4609  ELSE IF(metsol == 2) THEN
4610  WRITE(lu,*) ' METSOL = 2: diagonalization'
4611  ELSE IF(metsol == 3) THEN
4612  !GF WRITE(LU,*) ' METSOL = 3: MINRES'
4613  WRITE(lu,*) ' METSOL = 3: MINRES (rtol', mrestl,')'
4614  ELSE IF(metsol == 4) THEN
4615  WRITE(lu,*) ' METSOL = 4: GMRES'
4616  END IF
4617  WRITE(lu,*) ' with',mitera,' iterations'
4618  IF(matsto == 1) THEN
4619  WRITE(lu,*) ' MATSTO = 1: symmetric matrix, ', '(n*n+n)/2 elements'
4620  ELSE IF(matsto == 2) THEN
4621  WRITE(lu,*) ' MATSTO = 2: sparse matrix'
4622  END IF
4623  ! END IF
4624  IF(dflim /= 0.0) THEN
4625  WRITE(lu,103) 'Convergence assumed, if expected dF <',dflim
4626  END IF
4627 
4628  END DO ! print loop
4629 
4630  ! Wolfe conditions
4631 
4632  IF(0.0 < wolfc1.AND.wolfc1 < wolfc2.AND.wolfc2 < 1.0) go to 32
4633  IF(wolfc1 == 0.0) wolfc1=1.0e-4
4634  IF(wolfc2 == 0.0) wolfc2=0.9
4635  IF(0.0 < wolfc1.AND.wolfc1 < wolfc2.AND.wolfc2 < 1.0) go to 32
4636  IF(wolfc1 <= 0.0) wolfc1=1.0e-4
4637  IF(wolfc2 >= 1.0) wolfc2=0.9
4638  IF(wolfc1 > wolfc2) THEN ! exchange
4639  wolfc3=wolfc1
4640  wolfc1=wolfc2
4641  wolfc2=wolfc3
4642  ELSE
4643  wolfc1=1.0e-4
4644  wolfc2=0.9
4645  END IF
4646  WRITE(*,105) wolfc1,wolfc2
4647  WRITE(lun,105) wolfc1,wolfc2
4648 105 format(' Constants C1, C2 for Wolfe conditions:',g12.4,', ',g12.4)
4649 
4650  ! prepare matrix and gradient storage ------------------------------
4651  !32 CONTINUE
4652 32 matsiz(1)=int8(nagb)*int8(nagb+1)/2 ! number of words for double precision storage 'j'
4653  matsiz(2)=0 ! number of words for single precision storage '.'
4654  IF(matsto == 2) THEN ! sparse matrix
4655  matsiz(1)=ndimsa(3)+nagb
4656  matsiz(2)=ndimsa(4)
4657  CALL mpalloc(sparsematrixcolumns,ndimsa(2),'sparse matrix column list')
4658  CALL spbits(sparsematrixoffsets,sparsematrixcolumns,sparsematrixcompression)
4659  END IF
4660  matwords=matwords+matsiz(1)*2+matsiz(2) ! #words for matrix storage
4661 
4662  CALL feasma ! prepare constraint matrices
4663 
4664  CALL vmprep(matsiz) ! prepare matrix and gradient storage
4665  WRITE(*,*) ' '
4666  IF (matwords < 250000) THEN
4667  WRITE(*,*) 'Size of global matrix: < 1 MB'
4668  ELSE
4669  WRITE(*,*) 'Size of global matrix:',int(float(matwords)*4.0e-6),' MB'
4670  ENDIF
4671  ! print chi^2 cut tables
4672 
4673  ndfmax=naeqn-1
4674  WRITE(lunlog,*) ' '
4675  WRITE(lunlog,*) ' Cut values of Chi^2/Ndf and Chi2,'
4676  WRITE(lunlog,*) ' corresponding to 2 and 3 standard deviations'
4677  WRITE(lunlog,*) ' Ndf Chi^2/Ndf(2) Chi^2(2) ', &
4678  ' Chi^2/Ndf(3) Chi^2(3)'
4679  ndf=0
4680  DO
4681  IF(ndf > naeqn) exit
4682  IF(ndf < 10) THEN
4683  ndf=ndf+1
4684  ELSE IF(ndf < 20) THEN
4685  ndf=ndf+2
4686  ELSE IF(ndf < 100) THEN
4687  ndf=ndf+5
4688  ELSE IF(ndf < 200) THEN
4689  ndf=ndf+10
4690  ELSE
4691  exit
4692  END IF
4693  chin2=chindl(2,ndf)
4694  chin3=chindl(3,ndf)
4695  WRITE(lunlog,106) ndf,chin2,chin2*float(ndf),chin3, chin3*float(ndf)
4696  END DO
4697 
4698  WRITE(lunlog,*) 'LOOP2: ending'
4699  WRITE(lunlog,*) ' '
4700  CALL mend
4701 101 format(1x,a6,' =',i10,' = ',a)
4702 102 format(22x,a)
4703 103 format(1x,a,g12.4)
4704 106 format(i6,2(3x,f9.3,f12.1,3x))
4705 END SUBROUTINE loop2
4706 
4710 
4711 SUBROUTINE vmprep(msize)
4712  USE mpmod
4713  USE mpdalc
4714 
4715  IMPLICIT NONE
4716  INTEGER :: i
4717  INTEGER :: nmats
4718  !
4719  INTEGER(KIND=large), INTENT(IN) :: msize(2)
4720 
4721  INTEGER(KIND=large) :: length
4722  SAVE
4723  ! ...
4724  ! Vector/matrix storage
4725  length=nagb*mthrd
4726  CALL mpalloc(globalvector,length,'rhs vector') ! double precision vector
4727  lenglobalvec=nagb
4728  length=naeqn
4729  CALL mpalloc(localcorrections,length,'residual vector of one record')
4730  length=nalcn*nalcn
4731  CALL mpalloc(aux,length,' local fit scratch array: aux')
4732  CALL mpalloc(vbnd,length,' local fit scratch array: vbnd')
4733  CALL mpalloc(vbdr,length,' local fit scratch array: vbdr')
4734  length=((nalcn+1)*nalcn)/2
4735  CALL mpalloc(clmat,length,' local fit matrix: clmat')
4736  CALL mpalloc(vbk,length,' local fit scratch array: vbk')
4737  length=nalcn
4738  CALL mpalloc(blvec,length,' local fit vector: blvec')
4739  CALL mpalloc(vzru,length,' local fit scratch array: vzru')
4740  CALL mpalloc(scdiag,length,' local fit scratch array: scdiag')
4741  CALL mpalloc(scflag,length,' local fit scratch array: scflag')
4742  CALL mpalloc(ibandh,length,' local fit band width hist.: ibandh')
4743 
4744  CALL mpalloc(globalmatd,msize(1),'global matrix (D)' )
4745  CALL mpalloc(globalmatf,msize(2),'global matrix (F)')
4746 
4747  IF(metsol >= 3) THEN ! GMRES/MINRES algorithms
4748  ! array space is:
4749  ! variable-width band matrix or diagonal matrix for parameters
4750  ! followed by rectangular matrix for constraints
4751  ! followed by symmetric matrix for constraints
4752  IF(mbandw > 0) THEN ! variable-width band matrix
4753  length=nagb
4754  CALL mpalloc(indprecond,length,'pointer-array variable-band matrix')
4755  DO i=1,min(mbandw,nvgb)
4756  indprecond(i)=(i*i+i)/2 ! increasing number
4757  END DO
4758  DO i=min(mbandw,nvgb)+1,nvgb
4759  indprecond(i)=indprecond(i-1)+mbandw ! fixed band width
4760  END DO
4761  DO i=nvgb+1,nagb ! reset
4762  indprecond(i)=0
4763  END DO
4764  length=indprecond(nvgb)+ncgb*nvgb+(ncgb*ncgb+ncgb)/2
4765  CALL mpalloc(matprecond,length,'variable-band matrix')
4766  ELSE ! default preconditioner
4767  length=nvgb+ncgb*nvgb+(ncgb*ncgb+ncgb)/2
4768  CALL mpalloc(matprecond,length,'default preconditioner matrix')
4769  END IF
4770  END IF
4771 
4772 
4773  length=nagb
4774  CALL mpalloc(globalcorrections,length,'corrections') ! double prec corrections
4775 
4776  CALL mpalloc(workspacelinesearch,length,'auxiliary array (D2)') ! double aux 2
4777  CALL mpalloc(workspacei, length,'auxiliary array (I)') ! int aux 1
4778  CALL mpalloc(workspaceminres,length,'auxiliary array (D7)') ! double aux 7
4779 
4780  IF(metsol == 1) THEN
4781  CALL mpalloc(workspaced,length,'auxiliary array (D1)') ! double aux 1
4782  ! CALL MEGARR('t D',2*NAGB,'auxiliary array') ! double aux 8
4783  END IF
4784 
4785  IF(metsol == 2) THEN
4786  CALL mpalloc(workspaced,length,'auxiliary array (D1)') ! double aux 1
4787  CALL mpalloc(workspacediagonalization,length,'auxiliary array (D3)') ! double aux 3
4788  CALL mpalloc(workspaceeigenvalues,length,'auxiliary array (D6)') ! double aux 6
4789  length=nagb*nagb
4790  CALL mpalloc(workspaceeigenvectors,length,'(rotation) matrix U') ! rotation matrix
4791  END IF
4792 
4793  IF(metsol >= 3) THEN
4794  nmats=6
4795  length=nmats*nagb
4796  CALL mpalloc(workspaced,length,'auxiliary array (D1)') ! double aux 1
4797  END IF
4798 
4799 END SUBROUTINE vmprep
4800 
4804 
4805 SUBROUTINE minver
4806  USE mpmod
4807 
4808  IMPLICIT NONE
4809  INTEGER :: lun
4810  INTEGER :: ndefec
4811  INTEGER :: nrank
4812 
4813  SAVE
4814  ! ...
4815  lun=lunlog ! log file
4816  IF(lunlog == 0) lunlog=6
4817 
4818  ! WRITE(*,*) 'MINVER ICALCM=',ICALCM
4819  IF(icalcm == 1) THEN
4820  CALL sqminl(globalmatd, globalcorrections,nagb,nrank, &
4821  workspaced,workspacei)
4822  ndefec=nagb-nrank ! rank defect
4823  IF(ndefec /= 0) THEN
4824  WRITE(*,*) 'The rank defect of the symmetric',nagb, &
4825  '-by-',nagb,' matrix is ',ndefec,' (should be zero).'
4826  WRITE(lun,*) 'The rank defect of the symmetric',nagb, &
4827  '-by-',nagb,' matrix is ',ndefec,' (should be zero).'
4828  IF (iforce == 0) THEN
4829  isubit=1
4830  WRITE(*,*) ' --> enforcing SUBITO mode'
4831  WRITE(lun,*) ' --> enforcing SUBITO mode'
4832  END IF
4833  ELSE
4834  WRITE(lun,*) 'No rank defect of the symmetric matrix'
4835  END IF
4836 
4837  ELSE ! multiply gradient by inverse matrix
4838  CALL dbsvxl(globalmatd,globalvector,globalcorrections,nagb)
4839  END IF
4840 END SUBROUTINE minver
4841 
4843 SUBROUTINE mdiags
4844  USE mpmod
4845 
4846  IMPLICIT NONE
4847  REAL :: evalue
4848  INTEGER :: i
4849  INTEGER :: iast
4850  INTEGER :: idia
4851  INTEGER :: imin
4852  INTEGER :: lun
4853  INTEGER :: nmax
4854  INTEGER :: nmin
4855  INTEGER :: ntop
4856  INTEGER :: nvar
4857  !
4858  INTEGER(KIND=large) :: ii
4859  SAVE
4860  ! ...
4861 
4862  lun=lunlog ! log file
4863  IF(lunlog == 0) lun=6
4864 
4865  IF(icalcm == 1) THEN
4866  nvar=nagb
4867  DO i=1,nvar ! used in FEASIB
4868  ii=i
4869  workspaced(i)=globalmatd((ii*ii+ii)/2) ! save diagonal elements
4870  END DO
4871 
4872  ! eigenvalues eigenvectors symm_input
4873  CALL devrot(nvar,workspaceeigenvalues,workspaceeigenvectors,globalmatd, &
4874  workspacediagonalization,workspacei)
4875 
4876  ! histogram of positive eigenvalues
4877 
4878  nmax=int(1.0+log10(sngl(workspaceeigenvalues(1)))) ! > log of largest eigenvalue
4879  imin=1
4880  DO i=nagb,1,-1
4881  IF(workspaceeigenvalues(i) > 0.0d0) THEN
4882  imin=i ! index of smallest pos. eigenvalue
4883  exit
4884  END IF
4885  END DO
4886  nmin=int(log10(sngl(workspaceeigenvalues(imin)))) ! log of smallest pos. eigenvalue
4887  ntop=nmin+6
4888  DO WHILE(ntop < nmax)
4889  ntop=ntop+3
4890  END DO
4891 
4892  CALL hmpdef(7,float(nmin),float(ntop), 'log10 of positive eigenvalues')
4893  DO idia=1,nagb
4894  IF(workspaceeigenvalues(idia) > 0.0d0) THEN ! positive
4895  evalue=log10(sngl(workspaceeigenvalues(idia)))
4896  CALL hmpent(7,evalue)
4897  END IF
4898  END DO
4899  IF(nhistp /= 0) CALL hmprnt(7)
4900  CALL hmpwrt(7)
4901 
4902  iast=max(1,imin-60)
4903  CALL gmpdef(3,2,'low-value end of eigenvalues')
4904  DO i=iast,nagb
4905  evalue=sngl(workspaceeigenvalues(i))
4906  CALL gmpxy(3,float(i),evalue)
4907  END DO
4908  IF(nhistp /= 0) CALL gmprnt(3)
4909  CALL gmpwrt(3)
4910 
4911  DO i=1,nvar
4912  workspacediagonalization(i)=0.0d0
4913  IF(workspaceeigenvalues(i) /= 0.0d0) THEN
4914  workspacediagonalization(i)=max(0.0d0,log10(abs(workspaceeigenvalues(i)))+3.0d0)
4915  IF(workspaceeigenvalues(i) < 0.0d0) workspacediagonalization(i)=-workspacediagonalization(i)
4916  END IF
4917  END DO
4918  WRITE(lun,*) ' '
4919  WRITE(lun,*) 'The first (largest) eigenvalues ...'
4920  WRITE(lun,102) (workspaceeigenvalues(i),i=1,min(20,nagb))
4921  WRITE(lun,*) ' '
4922  WRITE(lun,*) 'The last eigenvalues ... up to',nvgb
4923  WRITE(lun,102) (workspaceeigenvalues(i),i=max(1,nvgb-19),nvgb)
4924  WRITE(lun,*) ' '
4925  WRITE(lun,*) 'The eigenvalues from',nvgb+1,' to',nagb
4926  WRITE(lun,102) (workspaceeigenvalues(i),i=nvgb+1,nagb)
4927  WRITE(lun,*) ' '
4928  WRITE(lun,*) 'Log10 + 3 of ',nagb,' eigenvalues in decreasing', ' order'
4929  WRITE(lun,*) '(for Eigenvalue < 0.001 the value 0.0 is shown)'
4930  WRITE(lun,101) (workspacediagonalization(i),i=1,nagb)
4931  IF(workspacediagonalization(nvar) < 0) WRITE(lun,*) 'Negative values are ', &
4932  'printed for negative eigenvalues'
4933  CALL devsig(nagb,workspaceeigenvalues,workspaceeigenvectors,globalvector,workspacediagonalization)
4934  WRITE(lun,*) ' '
4935  WRITE(lun,*) nvgb,' significances: insignificant if ', &
4936  'compatible with N(0,1)'
4937  WRITE(lun,101) (workspacediagonalization(i),i=1,nvgb)
4938 
4939 
4940 101 format(10f7.1)
4941 102 format(5e14.6)
4942 
4943  END IF
4944 
4945  ! solution ---------------------------------------------------------
4946 
4947  ! eigenvalues eigenvectors
4948  CALL devsol(nvar,workspaceeigenvalues,workspaceeigenvectors,globalvector,globalcorrections,workspacediagonalization)
4949  return
4950 END SUBROUTINE mdiags
4951 
4953 SUBROUTINE zdiags
4954  USE mpmod
4955 
4956  IMPLICIT NONE
4957  ! eigenvalue eigenvectors cov.matrix
4958  CALL devinv(nagb,workspaceeigenvalues,workspaceeigenvectors,globalmatd) ! inv
4959 
4960 END SUBROUTINE zdiags
4961 
4967 
4968 SUBROUTINE mminrs
4969  USE mpmod
4970 
4971  IMPLICIT NONE
4972  INTEGER :: istop
4973  INTEGER :: itn
4974  INTEGER :: itnlim
4975  INTEGER :: lun
4976  INTEGER :: nout
4977  INTEGER :: nrkd
4978  INTEGER :: nrkd2
4979 
4980  DOUBLE PRECISION :: shift
4981  DOUBLE PRECISION :: rtol
4982  DOUBLE PRECISION :: anorm
4983  DOUBLE PRECISION :: acond
4984  DOUBLE PRECISION :: rnorm
4985  DOUBLE PRECISION :: ynorm
4986  LOGICAL :: checka
4987  EXTERNAL avprod, mvsolv, mcsolv
4988  SAVE
4989  ! ...
4990  lun=lunlog ! log file
4991  IF(lunlog == 0) lun=6
4992 
4993  nout=lun
4994  itnlim=2000 ! iteration limit
4995  shift =0.0d0 ! not used
4996  rtol = mrestl ! from steering
4997  checka=.false.
4998 
4999  ! WRITE(*,*) 'METSOL MATSTO',METSOL,MATSTO
5000 
5001  ! WRITE(*,*) 'Input vector'
5002  ! WRITE(*,*) (DQ(IGVEC/2+I),I=1,NAGB)
5003 
5004  IF(mbandw == 0) THEN ! default preconditioner
5005  IF(icalcm == 1) THEN
5006  ! WRITE(LUN,*) 'MMINRS: PRECON started'
5007  CALL precon(ncgb,nvgb,matprecond,matprecond, matprecond(1+nvgb), &
5008  matprecond(1+nvgb+ncgb*nvgb))
5009  ! WRITE(LUN,*) 'MMINRS: PRECON ended'
5010  END IF
5011 
5012  ! WRITE(*,*) 'NAGB MBANDW ',NAGB,MBANDW
5013  ! WRITE(*,*) IGVEC,IDUX1,ISOLV,IDUX7
5014 
5015  CALL minres(nagb,globalvector, &
5016  workspaced(1), workspaced(1+nagb), workspaced(1+2*nagb), &
5017  workspaced(1+3*nagb),workspaced(1+4*nagb),workspaced(1+5*nagb), &
5018  globalcorrections,workspaceminres, avprod, mcsolv, checka ,.true. , shift, &
5019  nout , itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm)
5020  ! WRITE(*,*) 'MINRES ended'
5021  ! + DQ(IDUX1/2+1),DQ(IDUX8/2+1),DQ(IDUX3/2+1),
5022  ! + DQ(IDUX4/2+1),DQ(IDUX5/2+1),DQ(IDUX6/2+1),
5023  ELSE IF(mbandw > 0) THEN ! band matrix preconditioner
5024  IF(icalcm == 1) THEN
5025  WRITE(lun,*) 'MMINRS: EQUDEC started'
5026  CALL equdec(nvgb,ncgb,matprecond,indprecond,nrkd,nrkd2)
5027  WRITE(lun,*) 'MMINRS: EQUDEC ended'
5028  END IF
5029  CALL minres(nagb,globalvector, &
5030  workspaced(1), workspaced(1+nagb), workspaced(1+2*nagb), &
5031  workspaced(1+3*nagb),workspaced(1+4*nagb),workspaced(1+5*nagb), &
5032  globalcorrections,workspaceminres, avprod, mvsolv, checka ,.true. , shift, &
5033  nout, itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm)
5034  ELSE
5035  ! IF(ICALCM.NE.0) WRITE(LUN,*) 'MMINRS: no preconditioning !!!'
5036  CALL minres(nagb,globalvector, &
5037  workspaced(1), workspaced(1+nagb), workspaced(1+2*nagb), &
5038  workspaced(1+3*nagb),workspaced(1+4*nagb),workspaced(1+5*nagb), &
5039  globalcorrections,workspaceminres, avprod, mvsolv, checka ,.false., shift, &
5040  nout, itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm)
5041  END IF
5042  iitera=itn
5043  istopa=istop
5044  mnrsit=mnrsit+itn
5045 
5046  IF (istopa == 0) print *, 'MINRES: istop=0, exact solution x=0.'
5047 ! WRITE(*,*) 'Solution vector'
5048 ! WRITE(*,*) (DQ(ISOLV/2+I),I=1,NAGB)
5049 
5050 END SUBROUTINE mminrs
5051 
5059 
5060 SUBROUTINE mcsolv(n,x,y) ! solve M*y = x
5061  USE mpmod
5062 
5063  IMPLICIT NONE
5064  INTEGER,INTENT(IN) :: n
5065  DOUBLE PRECISION, INTENT(IN) :: x(n)
5066  DOUBLE PRECISION, INTENT(OUT) :: y(n)
5067  SAVE
5068  ! ...
5069  CALL presol(ncgb,nvgb,matprecond,matprecond(1+nvgb),matprecond(1+nvgb+ncgb*nvgb),y,x)
5070 END SUBROUTINE mcsolv
5071 
5079 
5080 SUBROUTINE mvsolv(n,x,y) ! solve M*y = x
5081  USE mpmod
5082 
5083  IMPLICIT NONE
5084  REAL :: factr
5085  INTEGER :: i
5086  INTEGER :: icgb
5087  INTEGER :: itgbi
5088  INTEGER :: ivgb
5089  INTEGER :: label
5090  INTEGER :: last
5091  INTEGER :: inone
5092 
5093  INTEGER, INTENT(IN) :: n
5094  DOUBLE PRECISION, INTENT(IN) :: x(n)
5095  DOUBLE PRECISION, INTENT(OUT) :: y(n)
5096  SAVE
5097  ! ...
5098  DO i=1,n
5099  y(i)=x(i) ! copy to output vector
5100  END DO
5101 
5102  DO icgb=1,ncgb ! copy previous lambda values
5103  ! Y(NVGB+ICGB)=DQ(ISOLV/2+NVGB+ICGB)
5104  END DO
5105 
5106  i=0
5107  icgb=0
5108  last=-1
5109  DO WHILE(i < lenconstraints)
5110  i=i+1
5111  label=listconstraints(i)%label
5112  IF(last == 0.AND.label == (-1)) icgb=icgb+1 ! next constraint header
5113  IF(label > 0) THEN
5114  factr=listconstraints(i)%value
5115  itgbi=inone(label) ! -> ITGBI= index of parameter label
5116  ivgb =globalparlabelindex(2,itgbi) ! -> variable-parameter index
5117  ! - A^T * lambda
5118  IF(ivgb /= 0) THEN
5119  ! Y(IVGB)=Y(IVGB) - DQ(ISOLV/2+NVGB+ICGB)*FACTR
5120  END IF
5121  END IF
5122  last=label
5123  END DO
5124 
5125  ! CALL VABSLV(NVGB,DQ(INDV/2+1),MQ(INDU+1),Y) ! solve for Y
5126  ! CALL VABSLV(NAGB,DQ(INDV/2+1),MQ(INDU+1),Y) ! solve for Y
5127 
5128  CALL equslv(nvgb,ncgb,matprecond,indprecond,y)
5129 END SUBROUTINE mvsolv
5130 
5131 
5132 
5133 !***********************************************************************
5134 
5147 
5148 SUBROUTINE xloopn !
5149  USE mpmod
5150 
5151  IMPLICIT NONE
5152  REAL :: catio
5153  REAL :: concu2
5154  REAL :: concut
5155  REAL, DIMENSION(2) :: ta
5156  INTEGER :: i
5157  INTEGER :: iact
5158  INTEGER :: iagain
5159  INTEGER :: idx
5160  INTEGER :: info
5161  INTEGER :: itgbi
5162  INTEGER :: ivgbi
5163  INTEGER :: jcalcm
5164  INTEGER :: k
5165  INTEGER :: labelg
5166  INTEGER :: litera
5167  INTEGER :: lrej
5168  INTEGER :: lun
5169  INTEGER :: lunp
5170  INTEGER :: minf
5171  INTEGER :: mrati
5172  INTEGER :: nan
5173  INTEGER :: nfaci
5174  INTEGER :: nloop
5175  INTEGER :: nrati
5176  INTEGER :: nrej
5177  INTEGER :: nsol
5178  INTEGER :: inone
5179 
5180  DOUBLE PRECISION :: stp
5181  DOUBLE PRECISION :: dratio
5182  DOUBLE PRECISION :: dfacin
5183  DOUBLE PRECISION :: djrat
5184  DOUBLE PRECISION :: dwmean
5185  DOUBLE PRECISION :: db
5186  DOUBLE PRECISION :: db1
5187  DOUBLE PRECISION :: db2
5188  DOUBLE PRECISION :: dbdot
5189  LOGICAL :: warner
5190  SAVE
5191  ! ...
5192 
5193  ! Printout of algorithm for solution and important parameters ------
5194 
5195  lun=lunlog ! log file
5196  IF(lunlog == 0) lunlog=6
5197 
5198  DO lunp=6,lunlog,lunlog-6
5199  WRITE(lunp,*) ' '
5200  WRITE(lunp,*) 'Solution algorithm: '
5201  WRITE(lunp,121) '=================================================== '
5202 
5203  IF(metsol == 1) THEN
5204  WRITE(lunp,121) 'solution method:','matrix inversion'
5205  ELSE IF(metsol == 2) THEN
5206  WRITE(lunp,121) 'solution method:','diagonalization'
5207  ELSE
5208  IF(metsol == 3) THEN
5209  WRITE(lunp,121) 'solution method:', 'minres (Paige/Saunders)'
5210  ELSE IF(metsol == 4) THEN
5211  WRITE(lunp,121) 'solution method:', &
5212  'gmres (generalized minimzation of residuals)'
5213  END IF
5214  END IF
5215  WRITE(lunp,123) 'convergence limit at Delta F=',dflim
5216  WRITE(lunp,122) 'maximum number of iterations=',mitera
5217  matrit=min(matrit,mitera)
5218  IF(matrit > 1) THEN
5219  WRITE(lunp,122) 'matrix recalculation up to ',matrit, '. iteration'
5220  END IF
5221  IF(metsol >= 3) THEN
5222  IF(matsto == 1) THEN
5223  WRITE(lunp,121) 'matrix storage:','full'
5224  ELSE IF(matsto == 2) THEN
5225  WRITE(lunp,121) 'matrix storage:','sparse'
5226  END IF
5227  WRITE(lunp,122) 'pre-con band-width parameter=',mbandw
5228  IF(mbandw == 0) THEN
5229  WRITE(lunp,121) 'pre-conditioning:','default'
5230  ELSE IF(mbandw < 0) THEN
5231  WRITE(lunp,121) 'pre-conditioning:','none!'
5232  ELSE IF(mbandw > 0) THEN
5233  WRITE(lunp,121) 'pre-conditioning=','band-matrix'
5234  END IF
5235  END IF
5236  IF(regpre == 0.0d0.AND.npresg == 0) THEN
5237  WRITE(lunp,121) 'using pre-sigmas:','no'
5238  ELSE
5239  ! FIXME: NPRESG contains parameters that failed the 'entries' cut...
5240  WRITE(lunp,124) 'pre-sigmas defined for', &
5241  float(100*npresg)/float(nvgb),' % of variable parameters'
5242  WRITE(lunp,123) 'default pre-sigma=',regpre
5243  END IF
5244  IF(nregul == 0) THEN
5245  WRITE(lunp,121) 'regularization:','no'
5246  ELSE
5247  WRITE(lunp,121) 'regularization:','yes'
5248  WRITE(lunp,123) 'regularization factor=',regula
5249  END IF
5250 
5251  IF(chicut /= 0.0) THEN
5252  WRITE(lunp,121) 'Chi square cut equiv 3 st.dev applied'
5253  WRITE(lunp,123) '... in first iteration with factor',chicut
5254  WRITE(lunp,123) '... in second iteration with factor',chirem
5255  WRITE(lunp,121) ' (reduced by sqrt in next iterations)'
5256  END IF
5257 
5258  IF(lhuber /= 0) THEN
5259  WRITE(lunp,122) 'Down-weighting of outliers in', lhuber,' iterations'
5260  WRITE(lunp,123) 'Cut on downweight fraction',dwcut
5261  END IF
5262 
5263 
5264 121 format(1x,a40,3x,a)
5265 122 format(1x,a40,2x,i2,a)
5266 123 format(1x,a40,2x,e9.2)
5267 124 format(1x,a40,3x,f5.1,a)
5268  END DO
5269 
5270  ! initialization of iterations -------------------------------------
5271 
5272  iitera=0
5273  nloop =0
5274  nsol =0 ! counter for solutions
5275  info =0
5276  lsinfo=0
5277  stp =0.0d0
5278  stepl =sngl(stp)
5279  concut=1.0e-12 ! initial constraint accuracy
5280  concu2=1.0e-06 ! constraint accuracy
5281  icalcm=1 ! require matrix calculation
5282  iterat=0 ! iteration counter
5283  iterat=-1
5284  litera=-2
5285  nrej=0 ! reset number of rejects
5286  IF(metsol == 1) THEN
5287  wolfc2=0.5 ! not accurate
5288  minf=1
5289  ELSE IF(metsol == 2) THEN
5290  wolfc2=0.5 ! not acurate
5291  minf=2
5292  ELSE IF(metsol == 3) THEN
5293  wolfc2=0.1 ! accurate
5294  minf=3
5295  ELSE IF(metsol == 4) THEN
5296  wolfc2=0.1 ! accurate
5297  minf=3
5298  END IF
5299 
5300  ! check initial feasibility of constraint equations ----------------
5301 
5302  WRITE(*,*) ' '
5303  IF(nofeas == 0) THEN ! make parameter feasible
5304  WRITE(lunlog,*) 'Checking feasibility of parameters:'
5305  WRITE(*,*) 'Checking feasibility of parameters:'
5306  CALL feasib(concut,iact) ! check feasibility
5307  IF(iact /= 0) THEN ! done ...
5308  WRITE(*,102) concut
5309  WRITE(*,*) ' parameters are made feasible'
5310  WRITE(lunlog,102) concut
5311  WRITE(lunlog,*) ' parameters are made feasible'
5312  ELSE ! ... was OK
5313  WRITE(*,*) ' parameters are feasible (i.e. satisfy constraints)'
5314  WRITE(lunlog,*) ' parameters are feasible (i.e. satisfy constraints)'
5315  END IF
5316  concut=concu2 ! cut for constraint check
5317  END IF
5318  iact=1 ! set flag for new data loop
5319  nofeas=0 ! set check-feasibility flag
5320 
5321  WRITE(*,*) ' '
5322  WRITE(*,*)'Reading files and accumulating vectors/matrices ...'
5323  WRITE(*,*) ' '
5324 
5325  CALL etime(ta,rstart)
5326  iterat=-1
5327  litera= 0
5328  jcalcm=-1
5329  iagain= 0
5330 
5331  icalcm=1
5332 
5333  ! Block 1: data loop with vector (and matrix) calculation ----------
5334 
5335  DO
5336  IF(iterat >= 0) THEN
5337  lcalcm=jcalcm+3 ! mode (1..4) of last loop
5338  IF(jcalcm+1 /= 0) THEN
5339  IF(iterat == 0) THEN
5340  CALL ploopa(6) ! header
5341  CALL ploopb(6)
5342  CALL ploopa(lunlog) ! iteration line
5343  CALL ploopb(lunlog)
5344  iterat=1
5345  CALL gmpxyd(1,float(nloopn),REAL(fvalue),0.5,0.) ! fcn-value graph (no Delta)
5346  ELSE
5347  IF(iterat /= litera) THEN
5348  CALL ploopb(6)
5349  ! CALL PLOOPA(LUNLOG)
5350  CALL ploopb(lunlog)
5351  litera=iterat
5352  CALL gmpxyd(1,float(nloopn),REAL(fvalue),0.5,delfun) ! fcn-value (with expected)
5353  IF(metsol == 3) THEN ! extend to 4, i.e. GMRES?
5354  CALL gmpxy(2,float(iterat),float(iitera)) ! MINRES iterations
5355  END IF
5356  ELSE
5357  CALL ploopc(6) ! sub-iteration line
5358  CALL ploopc(lunlog)
5359  CALL gmpxyd(1,float(nloopn),REAL(fvalue),0.5,0.) ! fcn-value graph (no Delta)
5360  END IF
5361  END IF
5362  ELSE
5363  CALL ploopd(6) ! solution line
5364  CALL ploopd(lunlog)
5365  END IF
5366  CALL etime(ta,rstart)
5367  ! CHK
5368  IF (iabs(jcalcm) <= 1) THEN
5369  idx=jcalcm+4
5370  times(idx )=(times(idx )*times(idx+3)+deltim) /(times(idx+3)+1.0)
5371  times(idx+3)= times(idx+3)+1.0
5372  END IF
5373  END IF
5374  jcalcm=icalcm
5375 
5376  IF(icalcm >= 0) THEN ! ICALCM = +1 & 0
5377  CALL loopn ! data loop
5378  CALL addcst ! constraints
5379  lrej=nrej
5380  nrej=nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) ! total number of rejects
5381  IF(3*nrej > nrecal) THEN
5382  WRITE(*,*) ' '
5383  WRITE(*,*) 'Data rejected in previous loop: '
5384  WRITE(*,*) ' ', &
5385  nrejec(0), ' (rank deficit/NaN) ',nrejec(1),' (Ndf=0) ', &
5386  nrejec(2), ' (huge) ',nrejec(3),' (large)'
5387  WRITE(*,*) 'Too many rejects (>33.3%) - stop'
5388  stop
5389  END IF
5390  END IF
5391  ! Block 2: new iteration with calculation of solution --------------
5392 
5393  IF(abs(icalcm) == 1) THEN ! ICALCM = +1 & -1
5394  DO i=1,nagb
5395  globalcorrections(i)=globalvector(i) ! copy rhs
5396  END DO
5397  DO i=1,nvgb
5398  itgbi=globalparvartototal(i)
5399  workspacelinesearch(i)=globalparameter(itgbi) ! copy X for line search
5400  END DO
5401  iterat=iterat+1 ! increase iteration count
5402  IF(metsol == 1) THEN
5403  CALL minver ! inversion
5404  ELSE IF(metsol == 2) THEN
5405  CALL mdiags ! diagonalization
5406  ELSE IF(metsol == 3) THEN
5407  CALL mminrs ! MINRES
5408  ELSE IF(metsol == 4) THEN
5409  WRITE(*,*) '... reserved for GMRES (not yet!)'
5410  CALL mminrs ! GMRES not yet
5411  END IF
5412 
5413  ! check feasibility and evtl. make step vector feasible
5414 
5415  DO i=1,nvgb
5416  itgbi=globalparvartototal(i)
5417  globalparcopy(itgbi)=globalparameter(itgbi) ! save
5418  globalparameter(itgbi)=globalparameter(itgbi)+globalcorrections(i) ! update
5419  END DO
5420  CALL feasib(concut,iact) ! improve constraints
5421  concut=concu2 ! new cut for constraint check
5422  DO i=1,nvgb
5423  itgbi=globalparvartototal(i)
5424  globalcorrections(i)=globalparameter(itgbi)-globalparcopy(itgbi) ! feasible stp
5425  globalparameter(itgbi)=globalparcopy(itgbi) ! restore
5426  END DO
5427 
5428  ! initialize line search based on slopes and prepare next
5429 
5430  CALL ptldef(wolfc2, 10.0, minf,10)
5431  IF(metsol == 1) THEN
5432  wolfc2=0.5 ! not accurate
5433  minf=3
5434  ELSE IF(metsol == 2) THEN
5435  wolfc2=0.5 ! not acurate
5436  minf=3
5437  ELSE IF(metsol == 3) THEN
5438  wolfc2=0.1 ! accurate
5439  minf=4
5440  ELSE IF(metsol == 4) THEN
5441  wolfc2=0.1 ! accurate
5442  minf=4
5443  END IF
5444  db=dbdot(nvgb,globalcorrections,globalvector)
5445  db1=dbdot(nvgb,globalcorrections,globalcorrections)
5446  db2=dbdot(nvgb,globalvector,globalvector)
5447  delfun=sngl(db)
5448  angras=sngl(db/dsqrt(db1*db2))
5449 
5450  IF(db <= 0.0d0) THEN
5451  WRITE(*,*) 'Function not decreasing:',db
5452  IF(db <= -1.0d-3) THEN ! 100311, VB/CK: allow some margin for numerics
5453  iagain=iagain+1
5454  IF (iagain <= 1) THEN
5455  WRITE(*,*) '... again matrix calculation'
5456  icalcm=1
5457  cycle
5458  ELSE
5459  WRITE(*,*) '... aborting iterations'
5460  go to 90
5461  END IF
5462  ELSE
5463  WRITE(*,*) '... stopping iterations'
5464  iagain=0
5465  go to 90
5466  END IF
5467  ELSE
5468  iagain=0
5469  END IF
5470  icalcm=0 ! switch
5471  ENDIF
5472  ! Block 3: line searching ------------------------------------------
5473 
5474  IF(icalcm+2 == 0) exit
5475  CALL ptline(nagb,workspacelinesearch, & ! current parameter values
5476  flines, & ! chi^2 function value
5477  globalvector, & ! gradient
5478  globalcorrections, & ! step vector stp
5479  stp, & ! returned step factor
5480  info) ! returned information
5481  ! WRITE(*,*) 'PTLINE returns INFO, STP=',INFO, STP
5482  lsinfo=info
5483 
5484  stepl=sngl(stp)
5485  nan=0
5486  DO i=1,nvgb
5487  itgbi=globalparvartototal(i)
5488  IF ((.NOT.(workspacelinesearch(i) <= 0.0d0)).AND. &
5489  (.NOT.(workspacelinesearch(i) > 0.0d0))) nan=nan+1
5490  globalparameter(itgbi)=workspacelinesearch(i) ! current parameter values
5491  END DO
5492 
5493  IF (nan > 0) THEN
5494  WRITE(*,*) 'Result vector containes ', nan,' NaNs - stop'
5495  stop
5496  END IF
5497 
5498  ! subito exit, if required -----------------------------------------
5499 
5500  IF(isubit /= 0) THEN ! subito
5501  WRITE(*,*) 'Subito! Exit after first step.'
5502  go to 90
5503  END IF
5504 
5505  IF(info == 0) THEN
5506  WRITE(*,*) 'INFO=0 should not happen (line search input err)'
5507  IF (iagain <= 0) THEN
5508  icalcm=1
5509  cycle
5510  ENDIF
5511  END IF
5512  IF(info < 0) cycle
5513  ! Block 4: line search convergence ---------------------------------
5514 
5515  CALL ptlprt(lunlog)
5516  CALL feasib(concut,iact) ! check constraints
5517  IF(iact /= 0.OR.chicut > 1.0) THEN
5518  icalcm=-1
5519  IF(iterat < matrit) icalcm=+1
5520  cycle ! iterate
5521  END IF
5522  IF(delfun <= dflim) go to 90 ! convergence
5523  IF(iterat >= mitera) go to 90 ! ending
5524  icalcm=-1
5525  IF(iterat < matrit) icalcm=+1
5526  cycle ! next iteration
5527 
5528  ! Block 5: iteration ending ----------------------------------------
5529 
5530 90 icalcm=-2
5531  END DO
5532  IF(nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) /= 0) THEN
5533  WRITE(*,*) ' '
5534  WRITE(*,*) 'Data rejected in last loop: '
5535  WRITE(*,*) ' ', &
5536  nrejec(0), ' (rank deficit/NaN) ',nrejec(1),' (Ndf=0) ', &
5537  nrejec(2), ' (huge) ',nrejec(3),' (large)'
5538  END IF
5539 
5540  dwmean=sumndf/dfloat(ndfsum)
5541  dratio=fvalue/dwmean/dfloat(ndfsum-nagb)
5542  catio=sngl(dratio)
5543  IF(lhuber /= 0) THEN
5544  IF (lhuber <= 3) THEN
5545  catio=catio/0.9326 ! correction Huber downweighting
5546  ELSE
5547  catio=catio/0.8228 ! correction Cauchy downweighting
5548  END IF
5549  END IF
5550  mrati=nint(100.0*catio)
5551 
5552  DO lunp=6,lunlog,lunlog-6
5553  WRITE(lunp,*) ' '
5554  IF (nfilw <= 0) THEN
5555  WRITE(lunp,*) 'Sum(Chi^2)/Sum(Ndf) =',fvalue
5556  WRITE(lunp,*) ' / (',ndfsum,'-',nagb,')'
5557  WRITE(lunp,*) ' =',dratio
5558  ELSE
5559  WRITE(lunp,*) 'Sum(W*Chi^2)/Sum(Ndf)/<W> =',fvalue
5560  WRITE(lunp,*) ' / (',ndfsum,'-', nagb,')'
5561  WRITE(lunp,*) ' /',dwmean
5562  WRITE(lunp,*) ' =',dratio
5563  END IF
5564  WRITE(lunp,*) ' '
5565  IF(lhuber /= 0) WRITE(lunp,*) &
5566  ' with correction for down-weighting ',catio
5567  END DO
5568  nrej=nrejec(0)+nrejec(1)+nrejec(2)+nrejec(3) ! total number of rejects
5569 
5570  ! ... the end with exit code ???????????????????????????????????????
5571 
5572  ! WRITE(*,199) ! write exit code
5573  ! + '-----------------------------------------------------------'
5574  ! IF(ITEXIT.EQ.0) WRITE(*,199)
5575  ! + 'Exit code = 0: Convergence reached'
5576  ! IF(ITEXIT.EQ.1) WRITE(*,199)
5577  ! + 'Exit code = 1: No improvement in last iteration'
5578  ! IF(ITEXIT.EQ.2) WRITE(*,199)
5579  ! + 'Exit code = 2: Maximum number of iterations reached'
5580  ! IF(ITEXIT.EQ.3) WRITE(*,199)
5581  ! + 'Exit code = 3: Failure'
5582  ! WRITE(*,199)
5583  ! + '-----------------------------------------------------------'
5584  ! WRITE(*,199) ' '
5585 
5586 
5587  nrati=nint(10000.0*float(nrej)/float(nrecal))
5588  djrat=0.01d0*dfloat(nrati)
5589  nfaci=nint(100.0*sqrt(catio))
5590 
5591  dratio=0.01d0*dfloat(mrati)
5592  dfacin=0.01d0*dfloat(nfaci)
5593 
5594  warner=.false.
5595  IF(mrati < 90.OR.mrati > 110) warner=.true.
5596  IF(nrati > 1) warner=.true.
5597  IF(nmiss1 /= 0) warner=.true.
5598  IF(iagain /= 0) warner=.true.
5599 
5600  IF(warner) THEN
5601  WRITE(*,199) ' '
5602  WRITE(*,199) ' '
5603  WRITE(*,199) 'WarningWarningWarningWarningWarningWarningWarningWarningWar'
5604  WRITE(*,199) 'arningWarningWarningWarningWarningWarningWarningWarningWarn'
5605  WRITE(*,199) 'rningWarningWarningWarningWarningWarningWarningWarningWarni'
5606  WRITE(*,199) 'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
5607  WRITE(*,199) 'ingWarningWarningWarningWarningWarningWarningWarningWarning'
5608  WRITE(*,199) 'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
5609  WRITE(*,199) 'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
5610 
5611  IF(mrati < 90.OR.mrati > 110) THEN
5612  WRITE(*,199) ' '
5613  WRITE(*,*) ' Chi^2/Ndf = ',dratio, ' (should be close to 1)'
5614  WRITE(*,*) ' => multiply all input standard ', &
5615  'deviations by factor',dfacin
5616  END IF
5617 
5618  IF(nrati > 1) THEN
5619  WRITE(*,199) ' '
5620  WRITE(*,*) ' Fraction of rejects =',djrat,' %', &
5621  ' (should be far below 1 %)'
5622  WRITE(*,*) ' => please provide correct mille data'
5623  END IF
5624 
5625  IF(iagain /= 0) THEN
5626  WRITE(*,199) ' '
5627  WRITE(*,*) ' Matrix not positiv definite '// &
5628  '(function not decreasing)'
5629  WRITE(*,*) ' => please provide correct mille data'
5630  END IF
5631 
5632  IF(nmiss1 /= 0) THEN
5633  WRITE(*,199) ' '
5634  WRITE(*,*) ' Rank defect =',nmiss1, &
5635  ' for constraint equations, should be 0'
5636  WRITE(*,*) ' => please correct constraint definition'
5637  END IF
5638 
5639  WRITE(*,199) ' '
5640  WRITE(*,199) 'WarningWarningWarningWarningWarningWarningWarningWarningWar'
5641  WRITE(*,199) 'arningWarningWarningWarningWarningWarningWarningWarningWarn'
5642  WRITE(*,199) 'rningWarningWarningWarningWarningWarningWarningWarningWarni'
5643  WRITE(*,199) 'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
5644  WRITE(*,199) 'ingWarningWarningWarningWarningWarningWarningWarningWarning'
5645  WRITE(*,199) 'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
5646  WRITE(*,199) 'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
5647  WRITE(*,199) ' '
5648 
5649  ENDIF
5650 
5651  CALL mend ! modul ending
5652 
5653  ! ------------------------------------------------------------------
5654 
5655  IF(metsol == 1) THEN
5656 
5657  ELSE IF(metsol == 2) THEN
5658  CALL zdiags
5659  ELSE IF(metsol == 3) THEN
5660  ! errors and correlations from MINRES
5661  DO k=1,mnrsel
5662  labelg=lbmnrs(k)
5663  IF(labelg == 0) cycle
5664  itgbi=inone(labelg)
5665  ivgbi=0
5666  IF(itgbi /= 0) ivgbi=globalparlabelindex(2,itgbi)
5667  IF(ivgbi < 0) ivgbi=0
5668  IF(ivgbi == 0) cycle
5669  ! determine error and global correlation for parameter IVGBI
5670  CALL solglo(ivgbi)
5671  END DO
5672 
5673  ELSE IF(metsol == 4) THEN
5674 
5675  END IF
5676 
5677  CALL prtglo ! print result
5678 
5679 102 format(' Call FEASIB with cut=',g10.3)
5680  ! 103 FORMAT(1X,A,G12.4)
5681 199 format(7x,a)
5682 END SUBROUTINE xloopn ! standard solution
5683 
5697 
5698 SUBROUTINE filetc
5699  USE mpmod
5700  USE mpdalc
5701 
5702  IMPLICIT NONE
5703  INTEGER :: i
5704  INTEGER :: ia
5705  INTEGER :: iargc
5706  INTEGER :: ib
5707  INTEGER :: ie
5708  INTEGER :: ierrf
5709  INTEGER :: ieq
5710  INTEGER :: ioff
5711  INTEGER :: iopt
5712  INTEGER :: ios
5713  INTEGER :: iosum
5714  INTEGER :: it
5715  INTEGER :: k
5716  INTEGER :: mat
5717  INTEGER :: nab
5718  INTEGER :: nline
5719  INTEGER :: npat
5720  INTEGER :: ntext
5721  INTEGER :: nu
5722  INTEGER :: nuf
5723  INTEGER :: nums
5724  INTEGER :: nufile
5725  INTEGER :: lenfileinfo
5726  INTEGER :: lenfilenames
5727  INTEGER :: matint
5728  INTEGER, DIMENSION(:,:), ALLOCATABLE :: vecfileinfo
5729  INTEGER, DIMENSION(:,:), ALLOCATABLE :: temparray
5730  INTEGER(KIND=large) :: rows
5731  INTEGER(KIND=large) :: cols
5732  INTEGER(KIND=large) :: newcols
5733  INTEGER(KIND=large) :: length
5734 
5735  CHARACTER (LEN=1024) :: text
5736  CHARACTER (LEN=1024) :: fname
5737  CHARACTER (LEN=14) :: bite(3)
5738  CHARACTER (LEN=32) :: keystx
5739  DOUBLE PRECISION :: dnum(100)
5740  SAVE
5741  DATA bite/'C_binary','text ','Fortran_binary'/
5742  ! ...
5743  CALL mstart('FILETC/X')
5744 
5745  nuf=1 ! C binary is default
5746  DO i=1,8
5747  times(i)=0.0
5748  END DO
5749 
5750  ! read command line options ----------------------------------------
5751 
5752  filnam=' ' ! print command line options and find steering file
5753  DO i=1,iargc()
5754  IF(i == 1) THEN
5755  WRITE(*,*) ' '
5756  WRITE(*,*) 'Command line options: '
5757  WRITE(*,*) '--------------------- '
5758  END IF
5759  CALL getarg(i,text) ! get I.th text from command line
5760  CALL rltext(text,ia,ib,nab) ! return indices for non-blank area
5761  WRITE(*,101) i,text(1:nab) ! echo print
5762  IF(text(ia:ia) /= '-') THEN
5763  nu=nufile(text(ia:ib)) ! inquire on file existence
5764  IF(nu == 2) THEN ! existing text file
5765  IF(filnam /= ' ') THEN
5766  WRITE(*,*) 'Second text file in command line - stop'
5767  stop
5768  ELSE
5769  filnam=text
5770  END IF
5771  ELSE
5772  WRITE(*,*) 'Open error for file:',text(ia:ib),' - stop'
5773  stop
5774  END IF
5775  ELSE
5776  IF(index(text(ia:ib),'b') /= 0) THEN
5777  mdebug=3 ! debug flag
5778  WRITE(*,*) 'Debugging requested'
5779  END IF
5780  it=index(text(ia:ib),'t')
5781  IF(it /= 0) THEN
5782  ictest=1 ! internal test files
5783  ieq=index(text(ia+it:ib),'=')+it
5784  IF (it /= ieq) THEN
5785  IF (index(text(ia+ieq:ib),'SL0' ) /= 0) ictest=2
5786  IF (index(text(ia+ieq:ib),'SLE' ) /= 0) ictest=3
5787  IF (index(text(ia+ieq:ib),'BP' ) /= 0) ictest=4
5788  IF (index(text(ia+ieq:ib),'BRLF') /= 0) ictest=5
5789  IF (index(text(ia+ieq:ib),'BRLC') /= 0) ictest=6
5790  END IF
5791  END IF
5792  IF(index(text(ia:ib),'s') /= 0) isubit=1 ! like "subito"
5793  IF(index(text(ia:ib),'f') /= 0) iforce=1 ! like "force"
5794  END IF
5795  IF(i == iargc()) WRITE(*,*) '--------------------- '
5796  END DO
5797 
5798 
5799  ! create test files for option -t ----------------------------------
5800 
5801  IF(ictest >= 1) THEN
5802  WRITE(*,*) ' '
5803  IF (ictest == 1) THEN
5804  CALL mptest ! 'wire chamber'
5805  ELSE
5806  CALL mptst2(ictest-2) ! 'silicon tracker'
5807  END IF
5808  IF(filnam == ' ') filnam='mp2str.txt'
5809  WRITE(*,*) ' '
5810  END IF
5811 
5812  ! check default steering file with file-name "steerfile" -----------
5813 
5814  IF(filnam == ' ') THEN ! check default steering file
5815  text='steerfile'
5816  CALL rltext(text,ia,ib,nab) ! return indices for non-blank area
5817  nu=nufile(text(ia:ib)) ! inquire on file existence and type
5818  IF(nu > 0) THEN
5819  filnam=text
5820  ELSE
5821  stop 'in FILETC: no steering file. .'
5822  END IF
5823  END IF
5824 
5825 
5826  ! open, read steering file:
5827  ! end
5828  ! fortranfiles
5829  ! cfiles
5830 
5831 
5832  CALL rltext(filnam,ia,ib,nfnam) ! return indices for non-blank area
5833  WRITE(*,*) ' '
5834  WRITE(*,*) 'Listing of steering file: ',filnam(1:nfnam)
5835  WRITE(*,*) '-------------------------'
5836  OPEN(10,file=filnam(1:nfnam),iostat=ios)
5837  IF(ios /= 0) THEN
5838  WRITE(*,*) 'Open error for steering file - stop'
5839  stop
5840  END IF
5841  ifile =0
5842  nfiles=0
5843 
5844  lenfileinfo=2
5845  lenfilenames=0
5846  rows=6; cols=lenfileinfo
5847  CALL mpalloc(vecfileinfo,rows,cols,'file info from steering')
5848  nline=0
5849  DO
5850  READ(10,102,iostat=ierrf) text ! read steering file
5851  IF (ierrf < 0) exit ! eof
5852  CALL rltext(text,ia,ib,nab) ! return indices for non-blank area
5853  nline=nline+1
5854  IF(nline <= 50) THEN ! print up to 50 lines
5855  WRITE(*,101) nline,text(1:nab)
5856  IF(nline == 50) WRITE(*,*) ' ...'
5857  END IF
5858 
5859  CALL rltext(text,ia,ib,nab) ! test content 'end'
5860  IF(ib == ia+2) THEN
5861  mat=matint(text(ia:ib),'end',npat,ntext)
5862  IF(mat == 3) THEN
5863  text=' '
5864  CALL intext(text,nline)
5865  WRITE(*,*) ' end-statement after',nline,' text lines'
5866  exit
5867  END IF
5868  END IF
5869 
5870  keystx='fortranfiles'
5871  !GF MAT=MATINT(TEXT(IA:IB),KEYSTX,NPAT,NTEXT)
5872  !GF IF(MAT.GE.NTEXT-NTEXT/10) THEN
5873  IF (text(ia:ib) == keystx) THEN ! exact comparison by GF
5874  nuf=3
5875  ! WRITE(*,*) 'Fortran files'
5876  cycle
5877  END IF
5878 
5879  keystx='Cfiles'
5880  !GF MAT=MATINT(TEXT(IA:IB),KEYSTX,NPAT,NTEXT)
5881  !GF IF(MAT.GE.NTEXT-NTEXT/10) THEN
5882  IF (text(ia:ib) == keystx) THEN ! exact comparison by GF
5883  nuf=1
5884  ! WRITE(*,*) 'Cfiles'
5885  cycle
5886  END IF
5887 
5888  ! file names
5889  ! check for file options (' -- ')
5890  ie=ib
5891  iopt=index(text(ia:ib),' -- ')
5892  IF (iopt > 0) ie=iopt-1
5893 
5894  IF(nab == 0) cycle
5895  nu=nufile(text(ia:ie)) ! inquire on file existence
5896  IF(nu > 0) THEN ! existing file
5897  IF (nfiles == lenfileinfo) THEN ! increase length
5898  CALL mpalloc(temparray,rows,cols,'temp file info from steering')
5899  temparray=vecfileinfo
5900  CALL mpdealloc(vecfileinfo)
5901  lenfileinfo=lenfileinfo*2
5902  newcols=lenfileinfo
5903  CALL mpalloc(vecfileinfo,rows,newcols,'file info from steering')
5904  vecfileinfo(:,1:cols)=temparray(:,1:cols)
5905  CALL mpdealloc(temparray)
5906  cols=newcols
5907  ENDIF
5908  nfiles=nfiles+1 ! count number of files
5909  IF(nu == 1) nu=nuf !
5910  lenfilenames=lenfilenames+ie-ia+1 ! total length of file names
5911  vecfileinfo(1,nfiles)=nline ! line number
5912  vecfileinfo(2,nfiles)=nu ! cbinary =1, text =2, fbinary=3
5913  vecfileinfo(3,nfiles)=ia ! file name start
5914  vecfileinfo(4,nfiles)=ie ! file name end
5915  vecfileinfo(5,nfiles)=iopt ! option start
5916  vecfileinfo(6,nfiles)=ib ! option end
5917  ELSE
5918  ! WRITE(*,*) 'Open error for file ',TEXT(IA:IB)
5919  ! STOP
5920  END IF
5921  END DO
5922  rewind 10
5923  ! read again to fill dynamic arrays with file info
5924  length=nfiles
5925  CALL mpalloc(mfd,length,'file type')
5926  CALL mpalloc(nfd,length,'file line (in steering)')
5927  CALL mpalloc(lfd,length,'file name length')
5928  CALL mpalloc(ofd,length,'file option')
5929  length=lenfilenames
5930  CALL mpalloc(tfd,length,'file name')
5931  nline=0
5932  i=1
5933  ioff=0
5934  DO
5935  READ(10,102,iostat=ierrf) text ! read steering file
5936  IF (ierrf < 0) exit ! eof
5937  nline=nline+1
5938  IF (nline == vecfileinfo(1,i)) THEN
5939  nfd(i)=vecfileinfo(1,i)
5940  mfd(i)=vecfileinfo(2,i)
5941  ia=vecfileinfo(3,i)-1
5942  lfd(i)=vecfileinfo(4,i)-ia ! length file name
5943  forall(k=1:lfd(i)) tfd(ioff+k)=text(ia+k:ia+k)
5944  ! tfd(i)=text(vecFileInfo(3,i):vecFileInfo(4,i)) ! file name
5945  ioff=ioff+lfd(i)
5946  ofd(i)=1.0 ! option for file
5947  IF (vecfileinfo(5,i) > 0) THEN
5948  CALL ratext(text(vecfileinfo(5,i)+4:vecfileinfo(6,i)),nums,dnum) ! translate text to DP numbers
5949  IF (nums > 0) ofd(i)=sngl(dnum(1))
5950  END IF
5951  i=i+1
5952  IF (i > nfiles) exit
5953  ENDIF
5954  ENDDO
5955  CALL mpdealloc(vecfileinfo)
5956  rewind 10
5957  ! additional info for binary files
5958  length=nfiles; rows=2
5959  CALL mpalloc(ifd,length,'integrated record numbers (=offset)')
5960  CALL mpalloc(jfd,length,'number of accepted records')
5961  CALL mpalloc(kfd,rows,length,'number of records in file, file order')
5962  CALL mpalloc(dfd,length,'ndf sum')
5963  CALL mpalloc(xfd,length,'max. record size')
5964  CALL mpalloc(wfd,length,'file weight')
5965  CALL mpalloc(cfd,length,'chi2 sum')
5966  !
5967  WRITE(*,*) '-------------------------'
5968  WRITE(*,*) ' '
5969 
5970  ! print table of files ---------------------------------------------
5971 
5972  IF (mprint > 1) THEN
5973  WRITE(*,*) 'Table of files:'
5974  WRITE(*,*) '---------------'
5975  END IF
5976  WRITE(8,*) ' '
5977  WRITE(8,*) 'Text and data files:'
5978  ioff=0
5979  DO i=1,nfiles
5980  forall(k=1:lfd(i)) fname(k:k)=tfd(ioff+k)
5981  ! fname=tfd(i)(1:lfd(i))
5982  IF (mprint > -1) WRITE(*,103) i,bite(mfd(i)),fname(1:lfd(i))
5983  WRITE(8,103) i,bite(mfd(i)),fname(1:lfd(i))
5984  ioff=ioff+lfd(i)
5985  END DO
5986  IF (mprint > 1) THEN
5987  WRITE(*,*) '---------------'
5988  WRITE(*,*) ' '
5989  END IF
5990 
5991  ! open the binary (data) files on unit 11, 12, ...
5992 
5993  iosum=0
5994  nfilf=0
5995  nfilb=0
5996  nfilw=0
5997  ioff=0
5998  DO i=1,nfiles ! Fortran files
5999  IF(mfd(i) == 3) THEN
6000  forall(k=1:lfd(i)) fname(k:k)=tfd(ioff+k)
6001  ! fname=tfd(i)(1:lfd(i))
6002  WRITE(*,*) 'Opening Fortran file ',10+nfilf+1, ' ',fname(1:lfd(i))
6003  OPEN(10+nfilf+1,file=fname(1:lfd(i)),iostat=ios, form='UNFORMATTED')
6004  IF(ios /= 0) THEN
6005  WRITE(*,*) 'Open error for file ',fname(1:lfd(i))
6006  iosum=iosum+1
6007  ELSE
6008  nfilf=nfilf+1
6009  nfilb=nfilb+1
6010  wfd(nfilb)=ofd(i)
6011  END IF
6012  END IF
6013  ioff=ioff+lfd(i)
6014  END DO
6015 
6016  ! open the binary C files
6017 
6018  nfilc=-1
6019  ioff=0
6020  DO i=1,nfiles ! Cfiles
6021  IF(mfd(i) == 1) THEN
6022 #ifdef READ_C_FILES
6023  IF(nfilc < 0) CALL initc(nfiles) ! uncommented by GF
6024  IF(nfilc < 0) nfilc=0
6025  forall(k=1:lfd(i)) fname(k:k)=tfd(ioff+k)
6026  ! fname=tfd(i)(1:lfd(i))
6027  WRITE(*,*) 'Opening C file ',nfilc+1, ': ',fname(1:lfd(i)) ! by GF
6028  CALL openc(fname(1:lfd(i)),ios)
6029  IF(ios /= 0) THEN
6030  WRITE(*,*) 'Open error ',ios,' for file ',fname(1:lfd(i))
6031  iosum=iosum+1 ! typo fixed by GF
6032  ELSE
6033  nfilc=nfilc+1
6034  nfilb=nfilb+1
6035  wfd(nfilb)=ofd(i)
6036  END IF
6037 #else
6038  WRITE(*,*) 'Opening of C-files not supported.'
6039  ! GF add
6040  iosum=iosum+1
6041  ! GF add end
6042 #endif
6043  END IF
6044  ioff=ioff+lfd(i)
6045  END DO
6046 
6047  DO k=1,nfilb
6048  kfd(1,k)=0 ! reset record counters
6049  kfd(2,k)=k ! set file number
6050  ifd(k)=0 ! reset integrated record numbers
6051  xfd(k)=0 ! reset max record size
6052  END DO
6053 
6054  IF(iosum /= 0) THEN
6055  stop 'FILETC: open error '
6056  END IF
6057  IF(nfilb == 0) THEN
6058  stop 'FILETC: no binary files '
6059  END IF
6060  WRITE(*,*) nfilb,' binary files opened' ! corrected by GF
6061 101 format(i3,2x,a)
6062 102 format(a)
6063 103 format(i3,2x,a14,3x,a)
6064  ! CALL mend
6065  return
6066 END SUBROUTINE filetc
6067 
6118 
6119 SUBROUTINE filetx ! ---------------------------------------------------
6120  USE mpmod
6121 
6122  IMPLICIT NONE
6123  INTEGER :: i
6124  INTEGER :: ia
6125  INTEGER :: ib
6126  INTEGER :: ierrf
6127  INTEGER :: ioff
6128  INTEGER :: ios
6129  INTEGER :: iosum
6130  INTEGER :: k
6131  INTEGER :: mat
6132  INTEGER :: nab
6133  INTEGER :: nfiln
6134  INTEGER :: nline
6135  INTEGER :: nlinmx
6136  INTEGER :: npat
6137  INTEGER :: ntext
6138  INTEGER :: matint
6139 
6140  ! CALL MSTART('FILETX')
6141 
6142  CHARACTER (LEN=1024) :: text
6143  CHARACTER (LEN=1024) :: fname
6144 
6145  WRITE(*,*) ' '
6146  WRITE(*,*) 'Processing text files ...'
6147  WRITE(*,*) ' '
6148 
6149  iosum=0
6150  ioff=0
6151  DO i=0,nfiles
6152  IF(i == 0) THEN
6153  WRITE(*,*) 'File ',filnam(1:nfnam)
6154  nlinmx=100
6155  ELSE
6156  nlinmx=10
6157  ia=ioff
6158  ioff=ioff+lfd(i)
6159  IF(mfd(i) /= 2) cycle ! exclude binary files
6160  forall(k=1:lfd(i)) fname(k:k)=tfd(ia+k)
6161  WRITE(*,*) 'File ',fname(1:lfd(i))
6162  IF (mprint > 1) WRITE(*,*) ' '
6163  OPEN(10,file=fname(1:lfd(i)),iostat=ios,form='FORMATTED')
6164  IF(ios /= 0) THEN
6165  WRITE(*,*) 'Open error for file ',fname(1:lfd(i))
6166  iosum=iosum+1
6167  cycle
6168  END IF
6169  END IF
6170 
6171  nline=0
6172  nfiln=1
6173  ! read text file
6174  DO
6175  READ(10,102,iostat=ierrf) text
6176  IF (ierrf < 0) THEN
6177  text=' '
6178  CALL intext(text,nline)
6179  WRITE(*,*) ' end-of-file after',nline,' text lines'
6180  exit ! eof
6181  ENDIF
6182  nline=nline+1
6183  IF(nline <= nlinmx.AND.mprint > 1) THEN ! print first 10 lines of every text fiLE
6184  CALL rltext(text,ia,ib,nab)
6185  nab=max(1,nab)
6186  WRITE(*,101) nline,text(1:nab)
6187  IF(nline == nlinmx) WRITE(*,*) ' ...'
6188  END IF
6189 
6190  CALL rltext(text,ia,ib,nab) ! test content 'end'
6191  IF(ib == ia+2) THEN
6192  mat=matint(text(ia:ib),'end',npat,ntext)
6193  IF(mat == 3) THEN
6194  text=' '
6195  CALL intext(text,nline)
6196  WRITE(*,*) ' end-statement after',nline,' text lines'
6197  exit
6198  END IF
6199  END IF
6200 
6201  IF(i == 0) THEN ! first text file - exclude lines with file names
6202  IF(nfiln <= nfiles.AND.nline == nfd(nfiln)) THEN
6203  nfiln=nfiln+1
6204  text=' '
6205  ! WRITE(*,*) 'line is excluded ',TEXT(1:10)
6206  END IF
6207  END IF
6208  ! WRITE(*,*) TEXT(1:40),' < interprete text'
6209  CALL intext(text,nline) ! interprete text
6210  END DO
6211  WRITE(*,*) ' '
6212  rewind 10
6213  CLOSE(unit=10)
6214  END DO
6215 
6216  IF(iosum /= 0) THEN
6217  stop 'FILETX: open error(s) in text files '
6218  END IF
6219 
6220  WRITE(*,*) '... end of text file processing.'
6221  WRITE(*,*) ' '
6222 
6223  IF(lunkno /= 0) THEN
6224  WRITE(*,*) ' '
6225  WRITE(*,*) lunkno,' unknown keywords in steering files, ', &
6226  'or file non-existing,'
6227  WRITE(*,*) ' see above!'
6228  WRITE(*,*) '------------> stop'
6229  WRITE(*,*) ' '
6230  stop
6231  END IF
6232 
6233  ! check methods
6234 
6235  IF(metsol == 0) THEN ! if undefined
6236  IF(matsto == 0) THEN ! if undefined
6237  ! METSOL=1 ! default is matrix inversion
6238  ! MATSTO=1 ! default is symmetric matrix
6239  ELSE IF(matsto == 1) THEN ! if symmetric
6240  metsol=3 ! MINRES
6241  ELSE IF(matsto == 2) THEN ! if sparse
6242  metsol=3 ! MINRES
6243  END IF
6244  ELSE IF(metsol == 1) THEN ! if inversion
6245  matsto=1 !
6246  ELSE IF(metsol == 2) THEN ! if diagonalization
6247  matsto=1
6248  ELSE IF(metsol == 3) THEN ! if MINRES
6249  ! MATSTO=2 or 1
6250  ELSE IF(metsol == 4) THEN ! if GMRES
6251  ! MATSTO=2 or 1
6252  ELSE
6253  WRITE(*,*) 'MINRES forced with sparse matrix!'
6254  WRITE(*,*) ' '
6255  WRITE(*,*) 'MINRES forced with sparse matrix!'
6256  WRITE(*,*) ' '
6257  WRITE(*,*) 'MINRES forced with sparse matrix!'
6258  metsol=3 ! forced
6259  matsto=2 ! forced
6260  END IF
6261  IF(matsto > 2) THEN
6262  WRITE(*,*) 'MINRES forced with sparse matrix!'
6263  WRITE(*,*) ' '
6264  WRITE(*,*) 'MINRES forced with sparse matrix!'
6265  WRITE(*,*) ' '
6266  WRITE(*,*) 'MINRES forced with sparse matrix!'
6267  metsol=3 ! forced
6268  matsto=2 ! forced
6269  END IF
6270 
6271  ! print information about methods and matrix storage modes
6272 
6273  WRITE(*,*) ' '
6274  WRITE(*,*) 'Solution method and matrix-storage mode:'
6275  IF(metsol == 1) THEN
6276  WRITE(*,*) ' METSOL = 1: matrix inversion'
6277  ELSE IF(metsol == 2) THEN
6278  WRITE(*,*) ' METSOL = 2: diagonalization'
6279  ELSE IF(metsol == 3) THEN
6280  WRITE(*,*) ' METSOL = 3: MINRES'
6281  ELSE IF(metsol == 4) THEN
6282  WRITE(*,*) ' METSOL = 4: GMRES (-> MINRES)'
6283 
6284  END IF
6285 
6286  WRITE(*,*) ' with',mitera,' iterations'
6287 
6288  IF(matsto == 1) THEN
6289  WRITE(*,*) ' MATSTO = 1: symmetric matrix, ', '(n*n+n)/2 elements'
6290  ELSE IF(matsto == 2) THEN
6291  WRITE(*,*) ' MATSTO = 2: sparse matrix'
6292  END IF
6293  IF(mbandw /= 0) THEN
6294  WRITE(*,*) ' and band matrix, width',mbandw
6295  END IF
6296 
6297 
6298 
6299  IF(chicut /= 0.0) THEN
6300  WRITE(*,*) 'Chi square cut equiv 3 st.dev applied ...'
6301  WRITE(*,*) ' in first iteration with factor',chicut
6302  WRITE(*,*) ' in second iteration with factor',chirem
6303  WRITE(*,*) ' (reduced by sqrt in next iterations)'
6304  END IF
6305 
6306  IF(lhuber /= 0) THEN
6307  WRITE(*,*) ' Down-weighting of outliers in', lhuber,' iterations'
6308  WRITE(*,*) ' Cut on downweight fraction',dwcut
6309  END IF
6310 
6311  CALL mend
6312 
6313 101 format(i3,2x,a)
6314 102 format(a)
6315 END SUBROUTINE filetx
6316 
6326 
6327 INTEGER FUNCTION nufile(fname)
6328 
6329  IMPLICIT NONE
6330  INTEGER :: ios
6331  INTEGER :: l1
6332  INTEGER :: ll
6333  INTEGER :: nm
6334  INTEGER :: npat
6335  INTEGER :: ntext
6336  INTEGER :: nuprae
6337  INTEGER :: matint
6338 
6339  CHARACTER (LEN=*), INTENT(INOUT) :: fname
6340  LOGICAL :: ex
6341  SAVE
6342  ! ...
6343  nufile=0
6344  IF(fname(1:5) == 'rfio:') nuprae=1
6345  IF(fname(1:5) == 'dcap:') nuprae=2
6346  IF(fname(1:5) == 'root:') nuprae=3
6347  IF(nuprae == 0) THEN
6348  INQUIRE(file=fname,iostat=ios,exist=ex)
6349  IF(ios /= 0) nufile=-abs(ios)
6350  IF(ios /= 0) return
6351  ELSE IF(nuprae == 1) THEN ! rfio:
6352  ll=len(fname)
6353  fname=fname(6:ll)
6354  ex=.true.
6355  nufile=1
6356  return
6357  ELSE
6358  ex=.true. ! assume file existence
6359  END IF
6360  IF(ex) THEN
6361  nufile=1 ! binary
6362  ll=len(fname)
6363  l1=max(1,ll-3)
6364  nm=matint('xt',fname(l1:ll),npat,ntext)
6365  IF(nm == 2) nufile=2 ! text
6366  IF(nm < 2) THEN
6367  nm=matint('tx',fname(l1:ll),npat,ntext)
6368  IF(nm == 2) nufile=2 ! text
6369  END IF
6370  END IF
6371 END FUNCTION nufile
6372 
6380 SUBROUTINE intext(text,nline)
6381  USE mpmod
6382  USE mptext
6383 
6384  IMPLICIT NONE
6385  INTEGER :: i
6386  INTEGER :: ia
6387  INTEGER :: ib
6388  INTEGER :: icount
6389  INTEGER :: ier
6390  INTEGER :: iomp
6391  INTEGER :: k
6392  INTEGER :: kkey
6393  INTEGER :: label
6394  INTEGER :: lkey
6395  INTEGER :: mat
6396  INTEGER :: miter
6397  INTEGER :: mxcnt
6398  INTEGER :: nab
6399  INTEGER :: nkey
6400  INTEGER :: nkeys
6401  INTEGER :: nl
6402  INTEGER :: nmeth
6403  INTEGER :: npat
6404  INTEGER :: ntext
6405  INTEGER :: nums
6406  INTEGER :: matint
6407 
6408  CHARACTER (LEN=*), INTENT(IN) :: text
6409  INTEGER, INTENT(IN) :: nline
6410 
6411  parameter(nkeys=9,nmeth=4)
6412  CHARACTER (LEN=16) :: methxt(nmeth)
6413  CHARACTER (LEN=16) :: keylst(nkeys)
6414  CHARACTER (LEN=32) :: keywrd
6415  CHARACTER (LEN=32) :: keystx
6416  DOUBLE PRECISION :: dnum(100)
6417  REAL :: plvs(3) ! vector array: real and ...
6418  INTEGER :: lpvs(3) ! ... integer
6419  equivalence(plvs(1),lpvs(1))
6420 
6421  INTERFACE
6422  SUBROUTINE additem(length,list,label,value)
6423  USE mpmod
6424  INTEGER, INTENT(IN OUT) :: length
6425  TYPE(listitem), DIMENSION(:), INTENT(IN OUT), ALLOCATABLE :: list
6426  INTEGER, INTENT(IN) :: label
6427  REAL, INTENT(IN) :: value
6428  END SUBROUTINE additem
6429  END INTERFACE
6430 
6431  DATA keylst/'unknown','parameter','constraint','wconstraint', &
6432  'measurement', 'method', &
6433  'mestimate', 'atleast','option'/
6434 
6435  !add number of iterations
6436  ! wconstraint like constraint
6437  ! measured r sigma label factor ...(more)
6438 
6439  SAVE
6440  DATA methxt/'diagonalization','inversion','fullMINRES', 'sparseMINRES'/
6441  DATA lkey/-1/ ! last keyword
6442 
6443  ! ...
6444  nkey=-1 ! new keyword
6445  CALL rltext(text,ia,ib,nab) ! return indices for non-blank area
6446  IF(nab == 0) goto 10
6447  CALL ratext(text(1:nab),nums,dnum) ! translate text to DP numbers
6448 
6449  IF(nums /= 0) nkey=0
6450  IF(keyb /= 0) THEN
6451  keywrd=text(keya:keyb) ! text is TEXT(KEYA)...TEXT(KEYB)
6452  ! WRITE(*,*) 'Keyword is ',KEYWRD
6453 
6454  ! compare keywords
6455 
6456  DO nkey=2,nkeys-1 ! loop over all pede keywords
6457  keystx=keylst(nkey) ! copy NKEY.th pede keyword
6458  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6459  IF(mat >= ntext-ntext/5) go to 10
6460  END DO
6461 
6462  ! more comparisons
6463 
6464  keystx='print'
6465  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6466  ! WRITE(*,*) KEYSTX,MAT,NTEXT
6467  ! IF(MAT.GE.NTEXT) THEN
6468  IF(mat >= (npat-npat/5)) THEN
6469  ! IF(MAT.GE.(NTEXT+NTEXT+3)/3) THEN
6470  mprint=1
6471  IF(nums > 0) mprint=idnint(dnum(1))
6472  return
6473  END IF
6474 
6475  keystx='debug'
6476  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6477  IF(mat >= (npat-npat/5)) THEN
6478  ! IF(MAT.GE.(NTEXT+NTEXT+3)/3) THEN
6479  mdebug=3
6480  ! GF IF(NUMS.GT.0) MPRINT=DNUM(1)
6481  IF(nums > 0) mdebug=idnint(dnum(1))
6482  IF(nums > 1) mdebg2=idnint(dnum(2))
6483  return
6484  END IF
6485 
6486  keystx='entries'
6487  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6488  IF(mat >= (npat-npat/5)) THEN
6489  ! IF(MAT.GE.(NTEXT+NTEXT+3)/3) THEN
6490  mreqen=0
6491  IF(nums > 0) mreqen=idnint(dnum(1))
6492  IF (mreqen <= 0) mreqen=1
6493  return
6494  END IF
6495 
6496  keystx='printrecord'
6497  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6498  IF(mat >= (npat-npat/5)) THEN
6499  IF(nums > 0) nrecpr=idnint(dnum(1))
6500  IF(nums > 1) nrecp2=idnint(dnum(2))
6501  return
6502  END IF
6503 
6504  keystx='maxrecord'
6505  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6506  IF(mat >= (npat-npat/5)) THEN
6507  IF (nums > 0.AND.dnum(1) > 0.) mxrec=idnint(dnum(1))
6508  return
6509  END IF
6510 
6511  keystx='cache'
6512  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6513  IF(mat >= (npat-npat/5)) THEN
6514  IF (nums > 0.AND.dnum(1) >= 0.) ncache=idnint(dnum(1)) ! cache size, <0 keeps default
6515  IF (nums == 2.AND.dnum(2) > 0..AND.dnum(2) <= 1.0) & ! read cache fill level
6516  fcache(1)=sngl(dnum(2))
6517  IF (nums >= 4) THEN ! explicit cache splitting
6518  DO k=1,3
6519  fcache(k)=sngl(dnum(k+1))
6520  END DO
6521  END IF
6522  return
6523  END IF
6524 
6525  keystx='chisqcut'
6526  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6527  IF(mat >= (npat-npat/5)) THEN
6528  IF(nums == 0) THEN ! always 3-sigma cut
6529  chicut=1.0
6530  chirem=1.0
6531  ELSE
6532  chicut=sngl(dnum(1))
6533  IF(chicut < 1.0) chicut=-1.0
6534  IF(nums == 1) THEN
6535  chirem=1.0 ! 3-sigma cut, if not specified
6536  ELSE
6537  chirem=sngl(dnum(2))
6538  IF(chirem < 1.0) chirem=1.0
6539  IF(chicut >= 1.0) chirem=min(chirem,chicut)
6540  END IF
6541  END IF
6542  return
6543  END IF
6544 
6545  ! GF added:
6546  keystx='hugecut'
6547  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6548  IF(mat >= (npat-npat/5)) THEN
6549  IF(nums > 0) chhuge=sngl(dnum(1))
6550  IF(chhuge < 1.0) chhuge=1.0 ! at least (!!) 3-sigma
6551  return
6552  END IF
6553  ! GF added end
6554 
6555  keystx='localfit'
6556  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6557  IF(mat >= (npat-npat/5)) THEN
6558  IF(nums > 0) lfitnp=idnint(dnum(1))
6559  IF(nums > 1) lfitbb=idnint(dnum(2))
6560  return
6561  END IF
6562 
6563  keystx='regularization'
6564  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6565  IF(mat >= (npat-npat/5)) THEN
6566  nregul=1
6567  regula=sngl(dnum(1))
6568  IF(nums >= 2) regpre=sngl(dnum(2))
6569  return
6570  END IF
6571 
6572  keystx='regularisation'
6573  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6574  IF(mat >= (npat-npat/5)) THEN
6575  nregul=1
6576  regula=sngl(dnum(1))
6577  IF(nums >= 2) regpre=sngl(dnum(2))
6578  return
6579  END IF
6580 
6581  keystx='presigma'
6582  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6583  IF(mat >= (npat-npat/5)) THEN
6584  regpre=sngl(dnum(1))
6585  return
6586  END IF
6587 
6588  keystx='matiter'
6589  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6590  IF(mat >= (npat-npat/5)) THEN
6591  matrit=idnint(dnum(1))
6592  return
6593  END IF
6594 
6595  keystx='matmoni'
6596  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6597  IF(mat >= (npat-npat/5)) THEN
6598  matmon=-1
6599  IF (nums > 0.AND.dnum(1) > 0.) matmon=idnint(dnum(1))
6600  return
6601  END IF
6602 
6603  keystx='bandwidth'
6604  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6605  ! IF(MAT.GE.(NTEXT+NTEXT+3)/3) THEN
6606  IF(mat >= (npat-npat/5)) THEN
6607  IF(nums > 0) mbandw=idnint(dnum(1))
6608  IF(mbandw < 0) mbandw=-1
6609  return
6610  END IF
6611 
6612  ! KEYSTX='outlierrejection'
6613  ! MAT=MATINT(TEXT(KEYA:KEYB),KEYSTX,NPAT,NTEXT) ! comparison
6614  ! WRITE(*,*) KEYSTX,MAT,(NTEXT+NTEXT)/3
6615  ! IF(MAT.GE.(NTEXT+NTEXT+NTEXT-2)/3) THEN
6616  ! IF(MAT.GE.(NPAT-NPAT/5)) THEN
6617  ! CHDFRJ=DNUM(1)
6618  ! IF(CHDFRJ.LT.3.0) CHDFRJ=100.0
6619  ! RETURN
6620  ! END IF
6621 
6622  ! KEYSTX='outliersuppression'
6623  ! MAT=MATINT(TEXT(KEYA:KEYB),KEYSTX,NPAT,NTEXT) ! comparison
6624  ! WRITE(*,*) KEYSTX,MAT,(NTEXT+NTEXT)/3
6625  ! IF(MAT.GE.(NTEXT+NTEXT+NTEXT-2)/3) THEN
6626  ! IF(MAT.GE.(NPAT-NPAT/5)) THEN
6627  ! LHUBER=DNUM(1)
6628  ! IF(LHUBER.LE.2) LHUBER=2 ! at least 2 Huber iterations
6629  ! RETURN
6630  ! END IF
6631 
6632  keystx='outlierdownweighting'
6633  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6634  ! WRITE(*,*) KEYSTX,MAT,(NTEXT+NTEXT)/3
6635  ! IF(MAT.GE.(NTEXT+NTEXT+NTEXT-2)/3) THEN
6636  IF(mat >= (npat-npat/5)) THEN
6637  lhuber=idnint(dnum(1))
6638  IF(lhuber > 0.AND.lhuber <= 2) lhuber=2 ! at least 2 Huber iterations (if any)
6639  return
6640  END IF
6641 
6642  keystx='dwfractioncut'
6643  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6644  IF(mat >= (npat-npat/5)) THEN
6645  dwcut=sngl(dnum(1))
6646  IF(dwcut > 0.5) dwcut=0.5
6647  return
6648  END IF
6649 
6650  keystx='pullrange'
6651  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6652  IF(mat >= (npat-npat/5)) THEN
6653  prange=abs(sngl(dnum(1)))
6654  return
6655  END IF
6656 
6657  keystx='subito'
6658  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6659  IF(mat >= (npat-npat/5)) THEN
6660  isubit=1
6661  return
6662  END IF
6663 
6664  keystx='force'
6665  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6666  IF(mat >= (npat-npat/5)) THEN
6667  iforce=1
6668  return
6669  END IF
6670 
6671  keystx='memorydebug'
6672  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6673  IF(mat >= (npat-npat/5)) THEN
6674  memdbg=1
6675  IF (nums > 0.AND.dnum(1) > 0.0) memdbg=idnint(dnum(1))
6676  return
6677  END IF
6678 
6679  keystx='globalcorr'
6680  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6681  IF(mat >= (npat-npat/5)) THEN
6682  igcorr=1
6683  return
6684  END IF
6685 
6686  keystx='threads'
6687  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6688  IF(mat >= (npat-npat/5)) THEN
6689  iomp=0
6690  !$ IOMP=1
6691  !$ IF (IOMP.GT.0) THEN
6692  !$ IF (NUMS.GE.1.AND.DNUM(1).GT.0.) MTHRD =IDNINT(DNUM(1))
6693  !$ IF (NUMS.GE.2) THEN
6694  !$ MTHRDR=MTHRD
6695  !$ IF (DNUM(2).GT.0.) MTHRDR=IDNINT(DNUM(2))
6696  !$ ENDIF
6697  !$ ELSE
6698  WRITE(*,*) 'WARNING: multithreading not available'
6699  !$ ENDIF
6700  return
6701  END IF
6702 
6703  keystx='compress'
6704  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6705  IF(mat >= (npat-npat/5)) THEN
6706  mcmprs=1
6707  return
6708  END IF
6709 
6710  keystx='errlabels'
6711  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6712  IF(mat >= (npat-npat/5).AND.mnrsel < 100) THEN
6713  nl=min(nums,100-mnrsel)
6714  DO k=1,nl
6715  lbmnrs(mnrsel+k)=idnint(dnum(k))
6716  END DO
6717  mnrsel=mnrsel+nl
6718  return
6719  END IF
6720 
6721  keystx='pairentries'
6722  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6723  IF(mat >= (npat-npat/5)) THEN
6724  ! This option could be implemented to get rid of parameter pairs
6725  ! that have very few entries - to save matrix memory size.
6726  IF (nums > 0.AND.dnum(1) > 0.0) THEN
6727  mreqpe=idnint(dnum(1))
6728  IF (nums >= 2.AND.dnum(2) >= dnum(1)) mhispe=idnint(dnum(2))
6729  IF (nums >= 3.AND.dnum(3) > 0.0) msngpe=idnint(dnum(3))
6730  icount=max(msngpe+1,mhispe)
6731  icount=max(mreqpe,icount)
6732  numbit=1 ! number of bits needed to count up to ICOUNT
6733  mxcnt=2
6734  DO k=1,30
6735  IF (icount >= mxcnt) THEN
6736  numbit=numbit+1
6737  mxcnt=mxcnt*2
6738  END IF
6739  END DO
6740  END IF
6741  return
6742  END IF
6743 
6744  keystx='wolfe'
6745  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6746  IF(mat >= (npat-npat/5)) THEN
6747  wolfc1=sngl(dnum(1))
6748  wolfc2=sngl(dnum(2))
6749  return
6750  END IF
6751 
6752  ! GF added:
6753  ! convergence tolerance for minres:
6754  keystx='mrestol'
6755  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6756  IF(mat >= (npat-npat/5)) THEN
6757  IF(nums > 0) THEN
6758  IF (dnum(1) < 1.0d-10.OR.dnum(1) > 1.0d-04) THEN
6759  WRITE(*,*) 'ERROR: need 1.0D-10 <= MRESTL ', &
6760  '<= 1.0D-04, but get ', dnum(1)
6761  ELSE
6762  mrestl=idnint(dnum(1))
6763  END IF
6764  END IF
6765  return
6766  END IF
6767  ! GF added end
6768 
6769  keystx='nofeasiblestart'
6770  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6771  IF(mat >= (npat-npat/5)) THEN
6772  nofeas=1 ! do not make parameters feasible at start
6773  return
6774  END IF
6775 
6776  keystx='histprint'
6777  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6778  IF(mat >= ntext-ntext/10) THEN
6779  ! IF(MAT.GE.(NPAT-NPAT/5)) THEN
6780  nhistp=1 ! print histograms
6781  return
6782  END IF
6783 
6784  keystx='fortranfiles'
6785  mat=matint(text(ia:ib),keystx,npat,ntext) ! comparison
6786  IF(mat >= ntext-ntext/10) return
6787 
6788  keystx='Cfiles'
6789  mat=matint(text(ia:ib),keystx,npat,ntext) ! comparison
6790  IF(mat >= ntext-ntext/10) return
6791 
6792  keystx=keylst(1)
6793  nkey=1 ! unknown keyword
6794  IF(nums /= 0) nkey=0
6795 
6796  WRITE(*,*) ' '
6797  WRITE(*,*) '**************************************************'
6798  WRITE(*,*) ' '
6799  WRITE(*,*) 'Unknown keyword(s): ',text(1:min(nab,50))
6800  WRITE(*,*) ' '
6801  WRITE(*,*) '**************************************************'
6802  WRITE(*,*) ' '
6803  lunkno=lunkno+1
6804 
6805  END IF
6806  ! result: NKEY = -1 blank
6807  ! NKEY = 0 numerical data, no text keyword or unknown
6808  ! NKEY > 0 keyword NKEY from list, keyword = KEYSTX
6809 
6810 
6811  ! content/lastcontent
6812  ! -------------------
6813  ! blank -1
6814  ! data 0
6815  ! keyword
6816  ! unknown 1
6817  ! parameter 2
6818  ! constraint 3
6819  ! method 4
6820  ! atleast 5
6821  ! option 6
6822  ! cut 7
6823 
6824 10 IF(nkey > 0) THEN ! new keyword
6825  lkey=nkey
6826  IF(lkey == 2) THEN ! parameter
6827  IF(nums == 3) THEN
6828  lpvs(1)=idnint(dnum(1)) ! label
6829  plvs(2)=sngl(dnum(2)) ! start value
6830  plvs(3)=sngl(dnum(3)) ! pre-sigma
6831  IF(lpvs(1) /= 0) THEN
6832  ! CALL megvec('7',plvs,3,0) ! always 3 words
6833  CALL additem(lenparameters,listparameters,lpvs(1),plvs(2))
6834  CALL additem(lenpresigmas,listpresigmas,lpvs(1),plvs(3))
6835  ELSE
6836  WRITE(*,*) 'Line',nline,' error, label=',lpvs(1)
6837  END IF
6838  ELSE IF(nums /= 0) THEN
6839  kkey=1 ! switch to "unknown" ?
6840  WRITE(*,*) 'Wrong text in line',nline
6841  WRITE(*,*) 'Status: new parameter'
6842  WRITE(*,*) '> ',text(1:nab)
6843  END IF
6844  ELSE IF(lkey == 3.OR.lkey == 4) THEN ! constraint
6845  ! WRITE(*,*) 'Keyword is constraint!',NUMS,' numerical data'
6846  IF(nums >= 1.AND.nums <= 2) THEN ! start constraint
6847  ! WRITE(*,*) 'NUMs is ',NUMS
6848  lpvs(1)=0
6849  plvs(2)=sngl(dnum(1)) ! r = r.h.s. value
6850  ! CALL megvec('8',plvs,2,0)
6851  CALL additem(lenconstraints,listconstraints,lpvs(1),plvs(2))
6852  ! WRITE(*,*) 'LPVS PLVS ',LPVS,PLVS
6853  lpvs(1)=-1 ! constraint
6854  IF(lkey == 4) lpvs(1)=-2 ! wconstraint (weighted)
6855  plvs(2)=0.0
6856  ! WRITE(*,*) 'LPVS PLVS ',LPVS,PLVS
6857  IF(nums == 2) plvs(2)=sngl(dnum(2)) ! sigma
6858  ! CALL megvec('8',plvs,2,0)
6859  CALL additem(lenconstraints,listconstraints,lpvs(1),plvs(2))
6860  ! WRITE(*,*) 'LPVS PLVS ',LPVS,PLVS
6861  ELSE
6862  kkey=1 ! switch to "unknown"
6863  WRITE(*,*) 'Wrong text in line',nline
6864  WRITE(*,*) 'Status: new keyword (w)constraint'
6865  WRITE(*,*) '> ',text(1:nab)
6866  END IF
6867  ELSE IF(lkey == 5) THEN ! measurement
6868  IF(nums == 2) THEN ! start measurement
6869  lpvs(1)=0
6870  plvs(2)=sngl(dnum(1)) ! r = r.h.s. value
6871  ! CALL megvec('9',plvs,2,0)
6872  CALL additem(lenmeasurements,listmeasurements,lpvs(1),plvs(2))
6873  lpvs(1)=-1 ! constraint
6874  plvs(2)=sngl(dnum(2)) ! sigma
6875  ! CALL megvec('9',plvs,2,0)
6876  CALL additem(lenmeasurements,listmeasurements,lpvs(1),plvs(2))
6877  ELSE
6878  kkey=1 ! switch to "unknown"
6879  WRITE(*,*) 'Wrong text in line',nline
6880  WRITE(*,*) 'Status: new keyword measurement'
6881  WRITE(*,*) '> ',text(1:nab)
6882  END IF
6883 
6884  ELSE IF(lkey == 6) THEN ! method
6885  miter=mitera
6886  IF(nums >= 1) miter=idnint(dnum(1))
6887  IF(miter >= 1) mitera=miter
6888  dflim=sngl(dnum(2))
6889  lkey=0
6890  DO i=1,nmeth
6891  keystx=methxt(i)
6892  mat=matint(text(keya:keyb),keystx,npat,ntext) ! comparison
6893  IF(mat >= ntext-ntext/5) THEN
6894  IF(i == 1) THEN ! diagonalization
6895  metsol=2
6896  matsto=1
6897  ELSE IF(i == 2) THEN ! inversion
6898  metsol=1
6899  matsto=1
6900  ELSE IF(i == 3) THEN ! fullMINRES
6901  metsol=3
6902  matsto=1
6903  ! METSOL=4 !!!
6904  ELSE IF(i == 4) THEN ! sparseMINRES
6905  ! METSOL=4
6906  metsol=3
6907  matsto=2
6908 
6909  END IF
6910  END IF
6911  END DO
6912  END IF
6913  ELSE IF(nkey == 0) THEN ! data for continuation
6914  IF(lkey == 2) THEN ! parameter
6915  IF(nums >= 3) THEN ! store data from this line
6916  lpvs(1)=idnint(dnum(1)) ! label
6917  plvs(2)=sngl(dnum(2)) ! start value
6918  plvs(3)=sngl(dnum(3)) ! pre-sigma
6919  IF(lpvs(1) /= 0) THEN
6920  ! CALL megvec('7',plvs,3,0) ! always 3 words
6921  CALL additem(lenparameters,listparameters,lpvs(1),plvs(2))
6922  CALL additem(lenpresigmas,listpresigmas,lpvs(1),plvs(3))
6923  ELSE
6924  WRITE(*,*) 'Line',nline,' error, label=',lpvs(1)
6925  END IF
6926  ELSE IF(nums > 1.AND.nums < 3) THEN
6927  kkey=1 ! switch to "unknown" ?
6928  WRITE(*,*) 'Wrong text in line',nline
6929  WRITE(*,*) 'Status continuation parameter'
6930  WRITE(*,*) '> ',text(1:nab)
6931  END IF
6932 
6933  ELSE IF(lkey == 3.OR.lkey == 4) THEN ! (w)constraint
6934  ier=0
6935  DO i=1,nums,2
6936  label=idnint(dnum(i))
6937  IF(label <= 0) ier=1
6938  END DO
6939  IF(mod(nums,2) /= 0) ier=1 ! reject odd number
6940  IF(ier == 0) THEN
6941  DO i=1,nums,2
6942  lpvs(1)=idnint(dnum(i)) ! label
6943  plvs(2)=sngl(dnum(i+1)) ! factor
6944  ! CALL megvec('8',plvs,2,0)
6945  CALL additem(lenconstraints,listconstraints,lpvs(1),plvs(2))
6946  END DO
6947  ELSE
6948  kkey=0
6949  WRITE(*,*) 'Wrong text in line',nline
6950  WRITE(*,*) 'Status continuation (w)constraint'
6951  WRITE(*,*) '> ',text(1:nab)
6952  END IF
6953 
6954  ELSE IF(lkey == 5) THEN ! measurement
6955  ! WRITE(*,*) 'continuation < ',NUMS
6956  ier=0
6957  DO i=1,nums,2
6958  label=idnint(dnum(i))
6959  IF(label <= 0) ier=1
6960  END DO
6961  IF(mod(nums,2) /= 0) ier=1 ! reject odd number
6962  ! WRITE(*,*) 'IER NUMS ',IER,NUMS
6963  IF(ier == 0) THEN
6964  DO i=1,nums,2
6965  lpvs(1)=idnint(dnum(i)) ! label
6966  plvs(2)=sngl(dnum(i+1)) ! factor
6967  ! CALL megvec('9',plvs,2,0)
6968  CALL additem(lenmeasurements,listmeasurements,lpvs(1),plvs(2))
6969  END DO
6970  ELSE
6971  kkey=0
6972  WRITE(*,*) 'Wrong text in line',nline
6973  WRITE(*,*) 'Status continuation measurement'
6974  WRITE(*,*) '> ',text(1:nab)
6975  END IF
6976 
6977  END IF
6978  END IF
6979 END SUBROUTINE intext
6980 
6988 SUBROUTINE additem(length,list,label,value)
6989  USE mpdalc
6990 
6991  INTEGER, INTENT(IN OUT) :: length
6992  TYPE(listitem), DIMENSION(:), INTENT(IN OUT), ALLOCATABLE :: list
6993  INTEGER, INTENT(IN) :: label
6994  REAL, INTENT(IN) :: value
6995 
6996  INTEGER(kind=large) :: newsize
6997  INTEGER(kind=large) :: oldsize
6998  TYPE(listitem), DIMENSION(:), ALLOCATABLE :: templist
6999 
7000  IF (label > 0.AND.value == 0.) return ! skip zero for valid labels
7001  IF (length == 0 ) THEN ! initial list with size = 100
7002  newsize = 100
7003  CALL mpalloc(list,newsize,' list ')
7004  ENDIF
7005  oldsize=size(list,kind=large)
7006  IF (length >= oldsize) THEN ! increase sizeby 20% + 100
7007  newsize = oldsize + oldsize/5 + 100
7008  CALL mpalloc(templist,oldsize,' temp. list ')
7009  templist=list
7010  CALL mpdealloc(list)
7011  CALL mpalloc(list,newsize,' list ')
7012  list(1:oldsize)=templist(1:oldsize)
7013  CALL mpdealloc(templist)
7014  ENDIF
7015  ! add to end of list
7016  length=length+1
7017  list(length)%label=label
7018  list(length)%value=value
7019 
7020 END SUBROUTINE additem
7021 
7023 SUBROUTINE mstart(text)
7024  USE mpmod, ONLY: textl
7025 
7026  IMPLICIT NONE
7027  INTEGER :: i
7028  INTEGER :: ka
7029  INTEGER :: kb
7030  INTEGER :: l
7031  CHARACTER (LEN=*), INTENT(IN) :: text
7032  CHARACTER (LEN=16) :: textc
7033  SAVE
7034  ! ...
7035  DO i=1,74
7036  textl(i:i)='_'
7037  END DO
7038  l=len(text)
7039  ka=(74-l)/2
7040  kb=ka+l-1
7041  textl(ka:kb)=text(1:l)
7042  WRITE(*,*) ' '
7043  WRITE(*,*) textl
7044  WRITE(*,*) ' '
7045  textc=text(1:l)//'-end'
7046 
7047  DO i=1,74
7048  textl(i:i)='_'
7049  END DO
7050  l=l+4
7051  ka=(74-l)/2
7052  kb=ka+l-1
7053  textl(ka:kb)=textc(1:l)
7054  return
7055 END SUBROUTINE mstart
7056 
7058 SUBROUTINE mend
7059  USE mpmod, ONLY: textl
7060 
7061  IMPLICIT NONE
7062  WRITE(*,*) ' '
7063  WRITE(*,*) textl
7064  CALL petime
7065  WRITE(*,*) ' '
7066 END SUBROUTINE mend
7067 
7074 
7075 SUBROUTINE mvopen(lun,fname)
7076  IMPLICIT NONE
7077  INTEGER :: l
7078  INTEGER, INTENT(IN) :: lun
7079  CHARACTER (LEN=*), INTENT(IN) :: fname
7080  CHARACTER (LEN=33) :: nafile
7081  CHARACTER (LEN=33) :: nbfile
7082  LOGICAL :: ex
7083  SAVE
7084  ! ...
7085  l=len(fname)
7086  IF(l > 32) stop 'File name too long '
7087  nafile=fname
7088  nafile(l+1:l+1)='~'
7089 
7090  INQUIRE(file=nafile(1:l),exist=ex)
7091  IF(ex) THEN
7092  INQUIRE(file=nafile(1:l+1),exist=ex)
7093  IF(ex) THEN
7094  CALL system('rm '//nafile)
7095  END IF
7096  nbfile=nafile
7097  nafile(l+1:l+1)=' '
7098  CALL system('mv '//nafile//nbfile)
7099  END IF
7100  OPEN(unit=lun,file=fname)
7101 END SUBROUTINE mvopen
7102 
7106 
7107 SUBROUTINE petime
7108  IMPLICIT NONE
7109  REAL :: ta(2)
7110  REAL :: rst
7111  REAL :: delta
7112  REAL :: rstp
7113  REAL :: secnd1
7114  REAL :: secnd2
7115  INTEGER :: ncount
7116  INTEGER :: nhour1
7117  INTEGER :: minut1
7118  INTEGER :: nsecd1
7119  INTEGER :: nhour2
7120  INTEGER :: minut2
7121  INTEGER :: nsecd2
7122 
7123  SAVE
7124  DATA ncount/0/
7125  ! ...
7126  ncount=ncount+1
7127  CALL etime(ta,rst)
7128  IF(ncount > 1) THEN
7129  delta=rst
7130  nsecd1=int(delta) ! -> integer
7131  nhour1=nsecd1/3600
7132  minut1=nsecd1/60-60*nhour1
7133  secnd1=delta-60*(minut1+60*nhour1)
7134  delta=rst-rstp
7135  nsecd2=int(delta) ! -> integer
7136  nhour2=nsecd2/3600
7137  minut2=nsecd2/60-60*nhour2
7138  secnd2=delta-60*(minut2+60*nhour2)
7139  WRITE(*,101) nhour1,minut1,secnd1, nhour2,minut2,secnd2
7140  END IF
7141 
7142  rstp=rst
7143  return
7144 101 format(i4,' h',i3,' min',f5.1,' sec total',18x,'elapsed', &
7145  i4,' h',i3,' min',f5.1,' sec')
7146 END SUBROUTINE petime ! print
7147 
7148 ! ----- accurate summation ----(from mpnum) ---------------------------------
7149 
7153 
7154 SUBROUTINE addsum(add)
7155  USE mpmod
7156 
7157  IMPLICIT NONE
7158  DOUBLE PRECISION:: add
7159  INTEGER::nadd
7160  ! ...
7161  nadd=int(add) ! convert to integer
7162  accuratensum=accuratensum+nadd ! sum integer
7163  accuratedsum=accuratedsum+(add-dfloat(nadd)) ! sum remainder
7164  IF(accuratedsum > 16.0d0) THEN ! + - 16
7165  accuratedsum=accuratedsum-16.0d0
7166  accuratensum=accuratensum+16
7167  END IF
7168  IF(accuratensum > nexp20) THEN ! if > 2^20: + - 2^20
7169  accuratenexp=accuratenexp+1
7170  accuratensum=accuratensum-nexp20
7171  END IF
7172  return
7173 END SUBROUTINE addsum
7174 
7178 
7179 SUBROUTINE getsum(asum)
7180  USE mpmod
7181 
7182  IMPLICIT NONE
7183  DOUBLE PRECISION, INTENT(OUT) ::asum
7184  asum=(accuratedsum+dfloat(accuratensum))+dfloat(accuratenexp)*dfloat(nexp20)
7185  accuratedsum=0.0d0
7186  accuratensum=0
7187  accuratenexp=0
7188  return
7189 END SUBROUTINE getsum
7190 
Generated on Wed May 23 2012 16:53:05 for Millepede-II by   doxygen 1.8.1